发点编程的东西,非局部均值法(NLM)降噪
图像降噪专用,不解释。

算法描述:

t1.jpg

t2.jpg

超慢的方法,处理每个像素的算法复杂度与处理半径的关系为O(N^4)。

偷懒一下加快算法:

t3.jpg

代码(写成了线程函数模式,以便多线程处理,NOISE为降噪强度,LERPC为lerp系数):



#include "stdafx.h"

#include <cmath>

#define NLM_WINDOW_RADIUS    3
#define NLM_BLOCK_RADIUS    3

#define NLM_WINDOW_AREA        ((2 * NLM_WINDOW_RADIUS + 1) * (2 * NLM_WINDOW_RADIUS + 1))
#define INV_NLM_WINDOW_AREA    (1.0f / (float)NLM_WINDOW_AREA)

#define NLM_WEIGHT_THRESHOLD    0.10f
#define NLM_LERP_THRESHOLD        0.10f

#define BLOCKDIM_X        8
#define BLOCKDIM_Y        8

__inline unsigned long vecLen(unsigned long p1, unsigned long p2)
{
    unsigned char *a = (unsigned char *)&p1;
    unsigned char *b = (unsigned char *)&p2;

    return ((b[0] - a[0]) * (b[0] - a[0])
          + (b[1] - a[1]) * (b[1] - a[1])
          + (b[2] - a[2]) * (b[2] - a[2]));
}

/*
__inline float lerpf(float a, float b, float c)
{
    return a + (b - a) * c;
}
*/

__inline unsigned long lerpf(unsigned long a, unsigned long b, long c)
{
    return a + (b - a) * c / 65536;
}


__inline unsigned long tex2D(unsigned long *image, int p, int h, int x, int y)
{
    if(x < 0)
        x = 0;
    else if(x >= p)
        x = p - 1;
    if(y < 0)
        y = 0;
    else if (y >= h)
        y = h - 1;

    return image[y * p + x];
}

#define NOISE        0.25f
#define LERPC        0.15f

UINT AFX_CDECL CPUFastNLMThread(LPVOID param)
{
    unsigned long *p = (unsigned long *)param;

    unsigned long *image = (unsigned long *)p[0];
    unsigned long pitch = p[1];
    unsigned long pitch4 = p[1] / 4;
    unsigned long width = p[2];
    unsigned long height = p[3];
    unsigned long y_amount = p[4];
    unsigned long y_start = p[5];
    unsigned long *image_out = (unsigned long *)p[6];

    float Noise = 1.0f / (NOISE * NOISE);
    float lerpC = LERPC;

    long bx, by, tx, ty;

    for(by = y_start; by < y_start + y_amount; by += BLOCKDIM_Y)
    {
        for(bx = 0; bx < width; bx += BLOCKDIM_Y)
        {
            // float fWeights[BLOCKDIM_X * BLOCKDIM_Y];
            long fWeights[BLOCKDIM_X * BLOCKDIM_Y];

            for(ty = 0; ty < BLOCKDIM_Y; ty ++)
                for(tx = 0; tx < BLOCKDIM_X; tx ++)
                {
                    long x = bx + tx;
                    long y = by + ty;
                    long cx = bx + NLM_WINDOW_RADIUS;
                    long cy = by + NLM_WINDOW_RADIUS;

                    if((x < width) && (y < height))
                    {
                        //Find color distance from current texel to the center of NLM window
                        // float weight = 0;
                        long weight = 0;
                        for(int n = -NLM_BLOCK_RADIUS; n <= NLM_BLOCK_RADIUS; n ++)
                            for(int m = -NLM_BLOCK_RADIUS; m <= NLM_BLOCK_RADIUS; m ++)
                                weight += vecLen(tex2D(image, pitch4, height, cx + m, cy + n), tex2D(image, pitch4, height, x + m, y + n));

                        //Geometric distance from current texel to the center of NLM window
                        float dist = (float)((tx - NLM_WINDOW_RADIUS) * (tx - NLM_WINDOW_RADIUS) + (ty - NLM_WINDOW_RADIUS) * (ty - NLM_WINDOW_RADIUS));

                        //Derive final weight from color and geometric distance
                        // weight = exp(-(weight / 65536.0f * Noise + dist * INV_NLM_WINDOW_AREA));

                        //Write the result to memory
                        // fWeights[ty * BLOCKDIM_X + tx] = weight;
                        fWeights[ty * BLOCKDIM_X + tx] = exp(-(weight / 65536.0f * Noise + dist * INV_NLM_WINDOW_AREA)) * 65536;
                    }
                }

            for(ty = 0; ty < BLOCKDIM_Y; ty ++)
                for(tx = 0; tx < BLOCKDIM_X; tx ++)
                {
                    long x = bx + tx;
                    long y = by + ty;
                    if((x < width) && (y < height))
                    {
                        //Normalized counter for the NLM weight threshold
                        // float fCount = 0;
                        long fCount = 0;
                        //Total sum of pixel weights
                        // float sumWeights = 0;
                        long sumWeights = 0;
                        //Result accumulator
                        // float clr[3] = {0, 0, 0};
                        long clr[3] = {0, 0, 0};
                        
                        int idx = 0;
                        
                        //Cycle through NLM window, surrounding (x, y) texel
                        for(long i = -NLM_WINDOW_RADIUS; i <= NLM_WINDOW_RADIUS + 1; i ++)
                            for(long j = -NLM_WINDOW_RADIUS; j <= NLM_WINDOW_RADIUS + 1; j ++)
                            {
                                //Load precomputed weight
                                // float weightIJ = fWeights[idx ++];
                                long weightIJ = fWeights[idx ++];
                                
                                //Accumulate (x + j, y + i) texel color with computed weight
                                long clrIJ = tex2D(image, pitch4, height, x + j, y + i);
                                clr[0] += (clrIJ & 0xff) * weightIJ;
                                clr[1] += ((clrIJ >> 8) & 0xff) * weightIJ;
                                clr[2] += ((clrIJ >> 16) & 0xff) * weightIJ;
                                
                                //Sum of weights for color normalization to [0..1] range
                                sumWeights  += weightIJ;
                                
                                //Update weight counter, if NLM weight for current window texel
                                //exceeds the weight threshold
                                // fCount += (weightIJ > NLM_WEIGHT_THRESHOLD) ? INV_NLM_WINDOW_AREA : 0;
                                fCount += (weightIJ > 6553) ? 1337 : 0;
                            }

                        //Normalize result color by sum of weights
                        /*
                        sumWeights = 1.0f / sumWeights;
                        clr[0] *= sumWeights;
                        clr[1] *= sumWeights;
                        clr[2] *= sumWeights;
                        */

                        clr[0] /= sumWeights;
                        clr[1] /= sumWeights;
                        clr[2] /= sumWeights;

                        //Choose LERP quotent basing on how many texels
                        //within the NLM window exceeded the weight threshold
                        // float lerpQ = (fCount > NLM_LERP_THRESHOLD * 65536) ? lerpC : 1.0f - lerpC;
                        long lerpQ = (fCount > 6553) ? 13107 : 52429;

                        //Write final result to memory
                        unsigned long clr00 = image[y * pitch4 + x];
                        clr[0] = lerpf(clr[0], clr00 & 0xff, lerpQ);
                        clr[1] = lerpf(clr[1], (clr00 >> 8) & 0xff, lerpQ);
                        clr[2] = lerpf(clr[2], (clr00 >> 16) & 0xff, lerpQ);
                        image_out[y * pitch4 + x] = ((unsigned char)clr[0]) | ((unsigned char)clr[1] << 8) | ((unsigned char)clr[2] << 16);
                    }
                }
        }
    }

    return 0;
}


+1000  科创币    phpskycn   2010-06-18    加少了。。。
+500  科创币    phpskycn   2010-06-18    不解释
来自:计算机科学 / 软件综合
 
小俊 作者
11年3个月前
1楼
为了提高速度,所有运算放大65536倍后基于整型进行,避免过多浮点运算。

运算速度在破双核E7200上大概为每秒处理160万像素,悲剧了,1000多万像素的照片得处理10秒。。。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
thallium2012
11年3个月前
2楼
你很厉害哈 呵呵 难以企及的高度 飘过
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊作者
11年3个月前
3楼
降噪效果测试:

500D拍的室内照片暗部,ISO1600,0.125s,原图裁剪:

TestDenoise.jpg

处理结果(参数为noise=0.25, lerpc=0.15)

TestDenoise_CPU.JPG
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊作者
11年3个月前
4楼
原算法来自CUDA SDK。

这个程序的CUDA版本代码在GTX285上的速度大概为双核E7200的80-90倍。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
y2k
11年3个月前
5楼
感覺細節失了不少...
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊作者
11年3个月前
6楼
引用第5楼y2k于2010-06-18 00:26发表的  :
感覺細節失了不少...


哪里的细节丢失了?请圈出

局部颜色过渡不自然倒是有点
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
y2k
11年3个月前
7楼
Cutt.JPG
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊作者
11年3个月前
8楼
引用第7楼y2k于2010-06-18 00:33发表的  :
Cutt.JPG


所有的花纹都还在吧,只是对比度损失了一些,可以加一级处理纠回来
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
delete
11年3个月前
9楼
处理后不锐了 不处理的图还可以接受
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊作者
11年3个月前
10楼
引用第9楼cat13于2010-06-18 00:46发表的  :
处理后不锐了 不处理的图还可以接受


世界上有一样算法叫USM锐化:

2.jpg

不处理那是噪点满天飞了。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
花落一天
11年3个月前
11楼
眼神不好,没看到燥点少了多少.........倒是看到锐度降低了很多,整个图糊了........
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
warmonkey
11年3个月前
12楼
对焦不太好,另外请上100%截图
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊作者
11年3个月前
13楼
引用第12楼warmonkey于2010-06-18 01:22发表的  :
对焦不太好,另外请上100%截图


已经是原图大小的截图,可能屏幕分辨率不够大的话论坛给缩小了

给个小尺寸截图:

3.jpg
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊作者
11年3个月前
14楼
引用第11楼花落一天于2010-06-18 01:19发表的  :
眼神不好,没看到燥点少了多少.........倒是看到锐度降低了很多,整个图糊了........


可能你的屏幕分辨率太低了,论坛把图缩小了
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
87685302
11年3个月前
15楼
那窗帘可以制作个裙子 那黑色的桌子让人肃穆起敬
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
ltl
11年3个月前
16楼
CUDA及其有爱……可惜缺钱缺硬件……
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
novakon
11年3个月前
17楼
这个方法去除了高频噪声,而且效果相当可以,只是低频部分仍然一塌糊涂。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
hmy0011
11年3个月前
18楼
用易语言的飘过。。。不懂高深的英文编程。。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊作者
11年3个月前
19楼
引用第18楼hmy0011于2010-06-18 13:31发表的  :
用易语言的飘过。。。不懂高深的英文编程。。


高不高深不在于语言本身,而在于算法。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
phpskycn
11年3个月前
20楼
看图
GTX285不是人人都有的,不过在GT240之类主流产品上的速度应该还可以接受。
上面那个E7200的速度是单线程的?另外感觉躁点仍然存在,少了很多 截图0010.jpg
截图0011.jpg
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊作者
11年3个月前
21楼
引用第20楼phpskycn于2010-06-18 15:02发表的  :
看图
GTX285不是人人都有的,不过在GT240之类主流产品上的速度应该还可以接受。
上面那个E7200的速度是单线程的?另外感觉躁点仍然存在,少了很多


是多线程的,线程数量等于CPU核的数量。我发的代码也是线程函数。

并行计算时代已经到来,过去靠彪CPU频率来拼速度的时代一去不复返。
即使不用CUDA,也应该在编程时考虑并行化和多线程。

另外GT240之流已经不能算是主流产品了。Fermi系列的最低端GF108已经达到GT240的级别。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
phpskycn
11年3个月前
22楼
GF108还有若干个月。。。
额那两个线程如何同步?把图像裁两半分开处理?
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊作者
11年3个月前
23楼
引用第22楼phpskycn于2010-06-18 17:10发表的  :
GF108还有若干个月。。。
额那两个线程如何同步?把图像裁两半分开处理?


正解。。。
GF108,1个月之内必出,96个SP。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
startwood
11年3个月前
24楼
额,那么中间分开那地方便还得再处理遍。。。
TSMC产能有点小小的问题,还有我发现NV似乎还没完全清仓,GT200系列还要存在一段时间,不知道啥时候GF10X才能替代原来那些产品,价格降到正常点的价位。
GF108要是便宜的话买一张……i7配7300LE太那个了。。。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊作者
11年3个月前
25楼
引用第24楼startwood于2010-06-18 17:35发表的  :
额,那么中间分开那地方便还得再处理遍。。。
TSMC产能有点小小的问题,还有我发现NV似乎还没完全清仓,GT200系列还要存在一段时间,不知道啥时候GF10X才能替代原来那些产品,价格降到正常点的价位。
GF108要是便宜的话买一张……i7配7300LE太那个了。。。


不需要处理分界部分。这个是根据输出像素来划分线程的。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
hypermania
11年3个月前
26楼
话说LZ编程看什么书?
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
startwood
11年3个月前
27楼
我个人认为LZ属于那种已经不用看书多年的。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
hypermania
11年3个月前
28楼
电子书是人类进步的电梯。。。一般的书都有电子书版本的是不。。。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊作者
11年3个月前
29楼
引用第27楼startwood于2010-06-18 19:28发表的  :
我个人认为LZ属于那种已经不用看书多年的。


纸质书的确很久没看过了,现在再看会睡着的。
电子文档倒是经常看,不一定是书,多数是技术文档或者论文一类的。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
ms1175555
11年3个月前
30楼
自觉原图记录了更多有用数据
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
wangzhm2020
11年1个月前
31楼
噪声少了,细节也损失了一些,有得必有失。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
boldness123
11年0个月前
32楼
马克?标记一下 免得找不找了
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
lovehongkong
10年10个月前
33楼
话说有没有RAW出jpeg时候直接降噪的??
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
cjc105808553
10年10个月前
34楼
楼主这么强大搞点计算机识别啊,生化危机中的计算机视觉很羡慕啊
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
phpskycn
10年8个月前
35楼
恕我挖坟...
今天尝试了N个降噪的软件,发现基本就是模糊化,再用USM后细节严重损失.
请问各位,有用NLM的现成的吗?
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
cwood
10年8个月前
36楼
回 19楼(小俊) 的帖子
其实处理的不错!挺好![s:251]
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
kaoheyton
8年3个月前
37楼
楼主能不能把代码发一份给我,最近在做并行图像处理,不知如何下手,想找个例子参考一下。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
evangil
8年3个月前
38楼
如楼上所说,细节丢失太多。
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
玄明
8年3个月前
39楼
回 38楼(kaoheyton) 的帖子
这样回复和发站内消息LZ都可以收到网站通知,悄悄的回个帖子谁也不知道你发帖了。论坛不欢迎以灌水的形式来挖坟,也不欢迎以个人原因来挖坟,请注意不要再犯咯
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
edowa
8年2个月前
40楼
这种效果很难接受,锐利度丢失严重,我宁愿噪点多点,看看Nokia Lumia 920的拍摄样张,它那个降噪算法才是真强大啊
回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
20!Dopaminor
3年0个月前
41楼

NVIDIA出的新算法

attachment icon Noise2Noise: Learning Image Restoration without Clean Data.pdf 9.45MB PDF 56次下载 预览


回复
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论

想参与大家的讨论?现在就 登录 或者 注册

所属专业
上级专业
同级专业
小俊
进士 学者 机友 笔友
文章
71
回复
1158
学术分
47
2006/12/29注册,8 天前活动
暂无简介
%7B%22isDisplay%22%3Atrue%7D

仅供内部学术交流或培训使用,请先保存到本地。本内容不代表科创观点,未经原作者同意,请勿转载。

插入资源
全部
图片
视频
音频
附件
全部
未使用
已使用
正在上传
空空如也~
上传中..{{f.progress}}%
处理中..
上传失败,点击重试
等待中...
{{f.name}}
空空如也~
(视频){{r.oname}}
{{selectedResourcesId.indexOf(r.rid) + 1}}
处理中..
处理失败
插入表情
我的表情
共享表情
Emoji
上传
注意事项
最大尺寸100px,超过会被压缩。为保证效果,建议上传前自行处理。
建议上传自己DIY的表情,严禁上传侵权内容。
点击重试等待上传{{s.progress}}%处理中...已上传
空空如也~
草稿箱
加载中...
此处只插入正文,如果要使用草稿中的其余内容,请点击继续创作。
{{fromNow(d.toc)}}
{{getDraftInfo(d)}}
标题:{{d.t}}
内容:{{d.c}}
继续创作
删除插入插入
{{forum.displayName}}
{{forum.countThreads}}
篇文章,
{{forum.countPosts}}
条回复
{{forum.description || "暂无简介"}}
ID: {{user.uid}}
学术分隐藏
{{submitted?"":"投诉或举报"}}
请选择违规类型:
{{reason.description}}
支持的图片格式:jpg, jpeg, png
插入公式
分享回复:{{shareId}}
加载中...
评论控制
加载中...
文号:{{pid}}
加载中...
详情
详情
推送到专栏从专栏移除
设为匿名取消匿名
查看作者
回复
只看作者
加入收藏取消收藏
加入关注取消关注
折叠回复
置顶取消置顶
评学术分
鼓励
设为精选取消精选
建议修改
编辑
通过审核
评论控制
退修或删除
历史版本
违规记录
投诉或举报
加入黑名单移除黑名单
查看IP
{{format('YYYY/MM/DD HH:mm:ss', toc)}}
下载资料
{{fileName}}
大小:{{size}}
下载当前附件将花费 {{costMessage}}
{{description}}
你当前剩余 {{holdMessage}}
{{fileName}}
大小:{{size}}
当前附件免费。
你已购买过此附件,下载当前附件不需要花费积分。
加载中...
{{errorInfo}}
附件已丢失
当前账号的附件下载数量限制如下:
时段 个数
{{f.startingTime}}点 - {{f.endTime}}点 {{f.fileCount}}