文档章节

O(1)效率的表面模糊算法优化。

abcijkxyz
 abcijkxyz
发布于 2016/11/22 16:41
字数 2200
阅读 4
收藏 0
点赞 0
评论 0

      本文的完整VS2013代码下载地址(解压密码本人博客名):http://files.cnblogs.com/files/Imageshop/SurfaceBlur.rar

      很久没有写文章了,主要是最近一段时间没有以前那么多空暇空间,内存和CPU占用率一致都很高,应前几日群里网友的要求,今天发个表面模糊的小程序来找回之前写博客的热情吧。

      国内我认为,破解表面模糊的原理的最早作者是我一直很崇拜的一位女士,她不会编程,英文也不怎么好,仅凭计算器和Excel两个工具破解了PS了很多  算法,真是个巾帼英雄。

      详见地址:http://www.missyuan.com/thread-428384-1-1.htm

      网上的有关该算法的matlab实现参考:http://www.cnblogs.com/tiandsp/archive/2012/11/06/2756441.html

      用C实现的参考:http://blog.csdn.net/maozefa/article/details/8270990

      表面模糊是属于典型的EPF滤波器中的一种,在PS的框架下好像也只有这一种自带的EPF算法,其核心也是数卷积的范畴,只是卷积的核是随着内容而变的,也属于方形半径内的算法,借助于直方图是可以做到于参数无关的O(1)算法。关于直方图的相关框架参考我的博文:任意半径局部直方图类算法在PC中快速实现的框架。, 但本文代码对其做了稍许改动。     

      为了表述方便,我们以灰度图像为例进行说明。首先,表面模糊有两个参数,半径Radius和阈值Threshold。 如果我们知道了以某点为中心,半径为Radius范围内的直方图数据Hist,以及该点的像素值,那根据原始的算法,其计算公式为:

//  最原始的算法
void Calc(unsigned short *Hist, unsigned char Value, int Threshold, unsigned char *&Pixel) { int Weight, Sum = 0, Divisor = 0; for (int Y = 0; Y < 256; Y++) { Weight = Hist[Y] * (2500 - abs(Y - Value) * 1000 / Threshold); if (Weight < 0) Weight = 0; Sum += Weight * Y; Divisor += Weight; } if (Divisor > 0) *Pixel = (Sum + (Divisor >> 1)) / Divisor; }

   注意这里我们为了减少浮点计算,将权重的计算公式放大了2500倍以便进行定点化,同时必须在最后增加一个Divisor > 0的判断,因为当Threshold很小时,可能会出现Divisor为0的现象。

     上述代码针对1000*1000的灰度图的执行时间约为1250ms,其中直方图的更新时间只有约50ms,速度难以接受。

     分析计算方法1,很明显权重计算的几个加减乘除以及下面的那句判断是比较耗时的,而其只是Y-Value的一个函数,因此,我们可以提前建立一个表,该表的索引范围从Min[Y - Value]到Max[Y - Value]之间,很明显,这个范围是[-255, 255],因此,建立如下的一个查找表:

 

for (int Y = -255; Y <= 255; Y++)
{
    int Factor = (2500 - abs(Y) * 1000 / Threshold);       
    if (Factor < 0) Factor = 0;
    Intensity[Y + 255] = Factor;
}

  有了这个查找表,我们来实现第二个版本的算法如下:

//    改进后的算法
unsigned char Calc2(unsigned short *Hist, unsigned char Value, unsigned short *Intensity)
{
    int Weight = 0, Sum = 0, Divisor = 0;
    unsigned short *Offset = Intensity + 255 - Value;
    for (int Y = 0; Y < 256; Y++)
    {
        Weight = Hist[Y] * Offset[Y];
        Sum += Weight * Y;
        Divisor += Weight;
    }
    if (Divisor > 0)
        return (Sum + (Divisor >> 1)) / Divisor;        //    四舍五入
    else
        return Value;
}

  同样大小的图,执行时间为350ms,速度提高约为3倍。

      我们接着来思考问题,上述有256个循环,如果我们将循环手动展开,会不会有提高呢, 我们先把代码更改如下:

//    优化后的算法
unsigned char Calc3(unsigned short *Hist, unsigned char Value, unsigned short *Intensity)       
{
    int Weight = 0, Sum = 0, Divisor = 0;
    unsigned short *Offset = Intensity + 255 - Value;
    Weight = Hist[0] * Offset[0];
    Sum += Weight * 0;  Divisor += Weight;        //    能不能用使用指令集的并行,没有去测试了
    Weight = Hist[1] * Offset[1];
    Sum += Weight * 1;  Divisor += Weight;
    Weight = Hist[2] * Offset[2];
    Sum += Weight * 2;  Divisor += Weight;
    Weight = Hist[3] * Offset[3];
    Sum += Weight * 3;  Divisor += Weight;
   /////////////////////////// ............................................................................
    Weight = Hist[251] * Offset[251];
    Sum += Weight * 251;  Divisor += Weight;
    Weight = Hist[252] * Offset[252];
    Sum += Weight * 252;  Divisor += Weight;
    Weight = Hist[253] * Offset[253];
    Sum += Weight * 253;  Divisor += Weight;
    Weight = Hist[254] * Offset[254];
    Sum += Weight * 254;  Divisor += Weight;
    Weight = Hist[255] * Offset[255];
    Sum += Weight * 255;  Divisor += Weight;
    if (Divisor > 0)
        return (Sum + (Divisor >> 1)) / Divisor;        //    四舍五入
    else
        return Value;
}

  为表述方便,中间省略了一些代码。

      测试结果为250ms,又快了一点点,为什么呢,我分析认为第一是减少了循环计数的时间,第二循环展开的 乘以 常数会被CPU优化为相关的移位或其他操作,而Calc2内部编译器是无法优化的。

      这样的函数系统一般是不会内联的,即使你在函数前面加上inline标识符,但是你可以在前面加上__forceinline标识,强制他内联,但是如果你这样做,你会发现速度反而会严重下降,为什么,请大家自行分析。

      我们在自己仔细看看,上面的循环很容易用SSE函数实现,既然我们的直方图的获取和更新利用了SSE,这里为什么不用呢,这样就诞生了我们的Calc4函数。

//    用SSE优化的算法
unsigned char Calc4(unsigned short *Hist, unsigned char Value, unsigned short *Intensity, unsigned short *Level)
{
    unsigned short *Offset = Intensity + 255 - Value;
    __m128i SumS = _mm_setzero_si128();
    __m128i WeightS = _mm_setzero_si128();
    for (int K = 0; K < 256; K += 8)
    {
        __m128i H = _mm_load_si128((__m128i const *)(Hist + K));
        __m128i L = _mm_load_si128((__m128i const *)(Level + K));                //    有能力可以使用256位的AVX寄存器
        __m128i I = _mm_loadu_si128((__m128i const *)(Offset + K));
        SumS = _mm_add_epi32(_mm_madd_epi16(_mm_mullo_epi16(L, I), H), SumS);
        WeightS = _mm_add_epi32(_mm_madd_epi16(H, I), WeightS);
    }
    const int *WW = (const int *)&WeightS;
    const int *SS = (const int *)&SumS;

    int Sum = SS[0] + SS[1] + SS[2] + SS[3];
    int Divisor = WW[0] + WW[1] + WW[2] + WW[3];
    if (Divisor > 0)
        return (Sum + (Divisor >> 1)) / Divisor;        //    四舍五入
    else
        return Value;
}

  关于上面几个SSE函数的使用,我不想多谈,也没啥难易理解的,注意其中的Level是我们为了方便,预定义的一个表,其形式如下:

for (int Y = 0; Y < 256; Y++)    Level[Y] = Y;            //    这个是为CalcSSE方便的使用的,其他两可以删除掉这里

     不定义这个也应该可以由其他的SSE函数构造k/k+1/k+2/k+3/k+4/k+5/k+6/k+7这样的__m128i变量,我这里这样做只是为了方便,你也可以自己更改下。

     我们直接把Calc4嵌入到程序中,运行,发现运行时间降低到了100ms,比Calc3有提高了2倍多,但是效果似乎不对,怎么回事呢。

     这主要是因为上述的SSE函数是针对unsigned short类型,而我们构造的Intensity数据较大,进行乘法后会超出unsigned short所能表达的范围,因此我们需要改动Intensity的定义:

    //    为了SSE里不溢出,把这里的数据变小,当然这样算法的准确度降低了,但是为了速度.......
    for (int Y = -255; Y <= 255; Y++)
    {
        int Factor = (255 - abs(Y) * 100 / Threshold);        
        if (Factor < 0) Factor = 0;
        Intensity[Y + 255] = Factor / 2;
    }

  最后一个除以2估计是因为SSE内部还是按照signed short处理的,这样做会导致算法的精度降低。

     经过上述改动,效果就正确了。

     对于彩色图像,一种做法就是直接扩展现在单通道的代码,让其支持三通道,另外一个办法就是把图像先拆分成3通道独立的数据,然后没通道独立处理,处理完成后再合成,这样做有两个好处,第一是代码复用;第二就是如果支持Openmp或者其他的并行库,可以让3通道并行起来执行。但是也有2个不足,第一是内存占用会增加很多,因为这种算法是不支持In-Place操作的,所以必须分配6份单通道的数据,而算法内部分配的内存由于并行的关系也要增加一些(不是三倍),及时考虑到可以把其中三个通道的数放置到Dest中,也会增加3份通道的数据,这对于某些设备可能是难以接受的(比如低端的安卓机)。具体如何使用就看应用场景了。 

      针对实际的应用,一种可选的进一步加速的方式就是把图像的色阶范围进一步缩小,比如由256色阶变为128或者64色阶,这样理论上还可以在快2倍到4倍,不过效果会稍有下降,一般128位时还是可以接受的。

 

         

     我看到很多人转载我的文章,我很感谢,但是很多人没有一点点的尊重别人的意识,转载请你在博文的最前面声明为转载,并不要更改本文下部赞助二维码。

    

 

    ****************************作者: laviewpbt   时间: 2015.10.24    联系QQ:  33184777 转载请保留本行信息**********************

本文转载自:http://www.cnblogs.com/Imageshop/p/5995093.html

共有 人打赏支持
abcijkxyz
粉丝 60
博文 6196
码字总数 1876
作品 0
深圳
项目经理
妹纸们的最爱 - 拓幻科技美颜SDK算法

现在各大手机制造商都在主推美颜效果,各种前后2000万像素,照亮你的美,各种逆光也清晰。其实这些看似神秘的美颜效果,除了依赖于手机像素之外,更重要的是攻城狮们对于美颜算法的构造。除了...

tillusory
01/31
0
0
高斯模糊算法的 C++ 实现

  2008 年在一个 PS 讨论群里,有网友不解 Photoshop 的高斯模糊中的半径是什么含义,因此当时我写了这篇文章:   对Photoshop高斯模糊滤镜的算法总结;   在那篇文章中,主要讲解了高...

hoodlum1980
2015/05/25
0
0
常用算法和复杂度总结

一、常用算法和复杂度 算法 名称 复杂度 备注 快速排序 QuickSort(A,p,r) 最坏:O(n2) 平均:O(nlog n) 均衡划分:O(nlog n) 合并排序 MergeSort(A,p,r) O(nlog n) 选最大 FindMax O(n) 选最...

啊莱
2010/01/03
0
0
RayTracking 光线跟踪算法

转自 http://www.cnblogs.com/daniagger/archive/2012/05/28/2521318.html 详细请看这位师兄的博客 1. Irradiance(辐照度) total amount of energy received per unit area of a surface 2.......

qq_41598072
04/13
0
0
【GPU精粹与Shader编程】(三) 《GPU Gems 1》全书核心内容提炼总结 · 下篇

本文由@浅墨_毛星云 出品,首发于知乎专栏,转载请注明出处 文章链接: https://zhuanlan.zhihu.com/p/36499291 题图背景来自《神秘海域4》。 系列文章前言 《GPU Gems》1~3 、《GPU Pro》1...

zhmxy555
05/06
0
0
干货 | 史上最全的支付宝二维码扫码优化技术方案

小蚂蚁说: 二维码又称二维条码,常见的二维码为QR Code,QR全称Quick Response,是一个近几年来移动设备上超流行的一种编码方式,它比传统的Bar Code条形码能存更多的信息,也能表示更多的数...

兔子酱
06/04
0
0
[发布] Photoshop 油画效果滤镜

    【原创性声明】本滤镜是我采用 PS SDK 开发而成,而滤镜的算法具体是由谁提出的可能不详,我是参考了 FilterExplorer 的源码(VC 6),本算法的主要参考来源是该项目中的 Filters.cp...

hoodlum1980
2011/01/15
0
0
哪些SQL语句会引起全表扫描

本文导读:大家都知道,用SQL语句对数据库进行操作时,如果引起全表扫描会对数据库的性能形成影响,下面向大家简单介绍SQL中哪些情况会引起全表扫描。 1、模糊查询效率很低: 原因:like本身...

zray4u
2016/10/18
81
0
时间复杂度和空间复杂度的理解

一个算法的时间复杂度其实就是这个算法跑过的次数 eg: 这个算法执行的总次数f(n)=n^2 + 2*n; 简单来说时间复杂度就是这个函数执行的基本操作次数 你是不是会想时间复杂度为什么不用时间来衡...

triorwy
2017/12/09
0
0
Java中的排序

以前总结过排序的种种知识 那么在Java中的Arrays.sort()是如何写的呢? JDK5中的Arrays.sort(int[]) JDK5基本类型的排序是使用优化了的快速排序,我们来看看JDK5中的优化点 /** * 将指定范围...

Hosee
2016/04/04
140
0

没有更多内容

加载失败,请刷新页面

加载更多

下一页

微信小程序Java登录流程(ssm实现具体功能和加解密隐私信息问题解决方案)

文章有不当之处,欢迎指正,如果喜欢微信阅读,你也可以关注我的微信公众号:好好学java,获取优质学习资源。 一、登录流程图 二、小程序客户端 doLogin:function(callback = () =>{}){let ...

公众号_好好学java
6分钟前
0
0
流利阅读笔记28-20180717待学习

“我不干了!” 英国脱欧大臣递交辞呈 雪梨 2018-07-17 1.今日导读 7 月 6 日,英国政府高官齐聚英国首相的官方乡间别墅——契克斯庄园,讨论起草了一份关于英国政府脱欧立场的白皮书。可是没...

aibinxiao
36分钟前
2
0
OSChina 周二乱弹 —— 理解超算排名这个事,竟然超出了很多人的智商

Osc乱弹歌单(2018)请戳(这里) 【今日歌曲】 @-冰冰棒- :分享Ed Sheeran/Beyoncé的单曲《Perfect Duet (with Beyoncé)》 《Perfect Duet (with Beyoncé)》- Ed Sheeran/Beyoncé 手机...

小小编辑
47分钟前
32
5
Android 获取各大音乐平台的真实下载地址

废话 电脑使用谷歌浏览器或者QQ浏览器的时候。。。。。。。说不清楚,还是看图吧 大概意思就是,只要网页上需要播放,只要能播放并且开始播放,这个过程就肯定会请求到相关的音乐资源,然后就...

她叫我小渝
今天
0
0
shell中的函数、shell中的数组、告警系统需求分析

shell中的函数 格式: 格式: function f_name() { command } 函数必须要放在最前面 示例1(用来打印参数) 示例2(用于定义加法) 示例3(用于显示IP) shell中的数组 shell中的数组1 定义数...

Zhouliang6
今天
2
0
用 Scikit-Learn 和 Pandas 学习线性回归

      对于想深入了解线性回归的童鞋,这里给出一个完整的例子,详细学完这个例子,对用scikit-learn来运行线性回归,评估模型不会有什么问题了。 1. 获取数据,定义问题     没有...

wangxuwei
今天
1
0
MAC安装MAVEN

一:下载maven压缩包(Zip或tar可选),解压压缩包 二:打开终端输入:vim ~/.bash_profile(如果找不到该文件新建一个:touch ./bash_profile) 三:输入i 四:输入maven环境变量配置 MAVEN_HO...

WALK_MAN
今天
0
0
33.iptables备份与恢复 firewalld的9个zone以及操作 service的操作

10.19 iptables规则备份和恢复 10.20 firewalld的9个zone 10.21 firewalld关于zone的操作 10.22 firewalld关于service的操作 10.19 iptables规则备份和恢复: ~1. 保存和备份iptables规则 ~2...

王鑫linux
今天
2
0
大数据教程(2.11):keeperalived+nginx高可用集群搭建教程

上一章节博主为大家介绍了目前大型互联网项目的系统架构体系,相信大家应该注意到其中很重要的一块知识nginx技术,在本节博主将为大家分享nginx的相关技术以及配置过程。 一、nginx相关概念 ...

em_aaron
今天
1
1
Apache Directory Studio连接Weblogic内置LDAP

OBIEE默认使用Weblogic内置LDAP管理用户及组。 要整理已存在的用户及组,此前办法是导出安全数据,文本编辑器打开认证文件,使用正则表达式获取用户及组的信息。 后来想到直接用Apache Dire...

wffger
今天
2
0

没有更多内容

加载失败,请刷新页面

加载更多

下一页

返回顶部
顶部