发点编程的东西,非局部均值法(NLM)降噪
小俊2010/06/18软件综合 IP:荷兰
图像降噪专用,不解释。

算法描述:

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 不解释
来自:计算机科学 / 软件综合
43
 
已屏蔽 原因:{{ notice.reason }}已屏蔽
{{notice.noticeContent}}
~~空空如也
小俊 作者
13年10个月前 IP:未同步
225468
为了提高速度,所有运算放大65536倍后基于整型进行,避免过多浮点运算。

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

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

TestDenoise.jpg

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

TestDenoise_CPU.jpg
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小俊作者
13年10个月前 IP:未同步
225473
原算法来自CUDA SDK。

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


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

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


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


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

2.jpg

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


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

给个小尺寸截图:

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


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


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


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

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

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


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


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


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

NVIDIA出的新算法

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


引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
amsl
2年4个月前 IP:山东
897782

这个算法不错效果还可以

如果可以一秒中连拍多张照片,剔除每张的早点,合成有效的像素,是不是就可以得到非常清晰的照片了.毕竟噪点是时刻变化的,每张的噪点也不在一个位置

引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
BSP
2年4个月前 IP:湖南
897799
引用amsl发表于42楼的内容
这个算法不错效果还可以如果可以一秒中连拍多张照片,剔除每张的早点,合成有效的像素,是不是就可以得到非...

但是每张照片之间的微小位移难以忽略,很难对齐完美

记得这个方法常用于固定摄影增强亮度,例如拍银河


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

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

所属专业
上级专业
同级专业
小俊
进士 学者 机友 笔友
文章
71
回复
1156
学术分
47
2006/12/29注册,16天18时前活动
暂无简介
主体类型:个人
所属领域:无
认证方式:手机号
IP归属地:未同步
文件下载
加载中...
{{errorInfo}}
{{downloadWarning}}
你在 {{downloadTime}} 下载过当前文件。
文件名称:{{resource.defaultFile.name}}
下载次数:{{resource.hits}}
上传用户:{{uploader.username}}
所需积分:{{costScores}},{{holdScores}}下载当前附件免费{{description}}
积分不足,去充值
文件已丢失

当前账号的附件下载数量限制如下:
时段 个数
{{f.startingTime}}点 - {{f.endTime}}点 {{f.fileCount}}
视频暂不能访问,请登录试试
仅供内部学术交流或培训使用,请先保存到本地。本内容不代表科创观点,未经原作者同意,请勿转载。
音频暂不能访问,请登录试试
支持的图片格式:jpg, jpeg, png
插入公式
评论控制
加载中...
文号:{{pid}}
投诉或举报
加载中...
{{tip}}
请选择违规类型:
{{reason.type}}

空空如也

加载中...
详情
详情
推送到专栏从专栏移除
设为匿名取消匿名
查看作者
回复
只看作者
加入收藏取消收藏
折叠回复
置顶取消置顶
评学术分
鼓励
设为精选取消精选
管理提醒
编辑
通过审核
评论控制
退修或删除
历史版本
违规记录
投诉或举报
加入黑名单移除黑名单
查看IP
{{format('YYYY/MM/DD HH:mm:ss', toc)}}