非局部均值降噪算法(NLM)原理及实现

news/2025/1/15 23:55:56/

文章目录

    • 一、概述
    • 二、算法原理
    • 三、算法流程
    • 四、MATLAB实现
    • 五、C++实现
    • 参考文献

一、概述

在日常生活中,最常见的 CT 图像噪声是高斯白噪声。目前,针对高斯白噪声的处理方法,主要有空间域中的以平滑为基本思想的均值滤波、高斯滤波、局部滤波等,此外还有频率域中的维纳滤波和小波阈值收缩等。局部强度统计特征是衡量区域内像素间的平均相似性,但这一特征难以准确辨别边缘与其邻近点之间的差异,导致了滤波结果中边缘信息的模糊。小波阈值收缩算法虽能够很好地估计信号的噪声,但是阈值和阈值函数的选择上存在不确定性,会造成信号的失真。
2005 年Buades 等提出了去除图像加性噪声的非局部均值滤波算法(Non-Local Means,NLM),该算法是利用图像非局部结构的相似性来去除噪声,恢复图像的主要几何结构。论文原文:A non-local algorithm for image denoising,还有一篇2011年的论文:Non-Local Means Denoising。

二、算法原理

在这里插入图片描述
如下所示为一般均值滤波的原理,其中权重仅通过中心像素和其邻域像素进行计算。当滤波半径过大的时候,图像的细节信息很难被保留,即使像双边滤波引入了像素间的灰度信息也会出现细节丢失严重的情况。
如下所示为双边滤波的公式及示意图:
在这里插入图片描述

非局部均值滤波是考虑到图像的非局部统计自相似性质,利用图像包含的大量重复结构进行去噪。通过计算非局部像素之间的相似度确定权值,然后加权平均来恢复图像。NLM 的主要特点 :将相似像素定义为具有相同邻域模式的像素,对像素周围固定大小的窗口内的信息进行比较,而不只是比较图像单个像素的信息,因此得到的相似性更加可靠。NLM基本原理如下 :
在这里插入图片描述
通过观察NLM的计算公式不难发现,其相对于常见的局部均值滤波,其区别主要体现在权重的计算上。在计算权重的相似度时,其通过像素间的邻域进行计算。

在这里插入图片描述
在这里插入图片描述
注意:通常非局部均值的滤波核都设置得比较大,如上图中所示的滤波核大小为9x9,因此NLM能够获得更佳的平滑度,同时更好的保留细节。

三、算法流程

在这里插入图片描述

四、MATLAB实现

无优化的实现:

function [output]=nlm_filter(input,t,f,h)% 输入: 待平滑的图像
% t: 搜索窗口半径
% f: 相似性窗口半径
% h: 平滑参数
% NLmeans(ima,5,2,sigma);% 图像大小
[m n]=size(input);
% 输出
output=zeros(m,n);
input2 = padarray(input,[f+t f+t],'symmetric');%边界作对称处理% 高斯核
ksize = f*2 + 1;
kernel = fspecial('gaussian',[ksize,ksize],1);
kernel = kernel / sum(sum(kernel));
kernel(:) = 1;
h=h*h;for i=1:mfor j=1:ni1 = i+ f+t;%原始图像的像素位置 (中心像素)j1 = j+ f+t;W1= input2(i1-f:i1+f , j1-f:j1+f);%小窗口wmax=0;average=0;sweight=0;%rmin = max(i1-t,f+1);%rmax = min(i1+t,m+f);%smin = max(j1-t,f+1);%smax = min(j1+t,n+f);rmin=i1-t;rmax=i1+t;smin=j1-t;smax=j1+t;for r=rmin:1:rmax %大窗口for s=smin:1:smaxif(r==i1 && s==j1)continue;endW2= input2(r-f:r+f , s-f:s+f);    %大搜索窗口中的小相似性窗口d = sum(sum(kernel.*(W1-W2).*(W1-W2)));w=exp(-d/h); %权重if w>wmaxwmax=w;   %求最大权重endsweight = sweight + w;  %大窗口中的权重和average = average + w*input2(r,s);endendaverage = average + wmax*input2(i1,j1);sweight = sweight + wmax;if sweight > 0output(i,j) = average / sweight;elseoutput(i,j) = input(i,j);endend
end

使用积分图加速(这部分还没仔细看其实现原理):

clc;
clear;
close all;src=imread('lena.jpg');
src=rgb2gray(src);
src=double(src);
figure,imshow(src,[]),title('src image')[m,n]=size(src);
ds=2;% block size for calculate weight
Ds=5;% search block
h=10;% decay factor
offset=ds+Ds;
PaddedImg = padarray(src,[ds+Ds,ds+Ds],'symmetric','both');% 扩展图像,便于处理sumimage=zeros(m,n);
sumweight=zeros(m,n);
maxweight=zeros(m,n);
image=PaddedImg(1+Ds:Ds+m+ds,1+Ds:Ds+n+ds);
[M,N]=size(image);
for r=-Ds:Dsfor s=-Ds:Ds%跳过当前点偏移if(r==0&&s==0)continue;end%求得差值积分图wimage=PaddedImg(1+Ds+r:Ds+m+ds+r,1+Ds+s:Ds+n+ds+s);diff=image-wimage;diff=diff.^2;J=cumsum(diff,1);J=cumsum(J,2);%计算距离distance=J(M-m+1:M,N-n+1:N)+J(1:m,1:n)-J(M-m+1:M,1:n)-J(1:m,N-n+1:N);distance=distance/((2*ds+1).^2);%计算权重并获得单个偏移下的加权图像weight=exp(-distance./(h*h));sumimage=sumimage+weight.*wimage(ds+1:ds+m,ds+1:ds+n);sumweight=sumweight+weight;maxweight=max(maxweight,weight);end
end
sumimage=sumimage+maxweight.*image(ds+1:ds+m,ds+1:ds+n);
sumweight=sumweight+maxweight;
dst=sumimage./sumweight;figure,imshow(dst,[]),title('dst');

降噪前后的效果比较:
在这里插入图片描述

五、C++实现

//计算0~255的平方查找表
float table1[256];
static void cal_lookup_table1(void)
{for (int i = 0; i < 256; i++){table1[i] = (float)(i*i);}
}//计算两个0~255的数的绝对差值的查找表
uchar table2[256][256];
static void cal_lookup_table2(void)
{for (int i = 0; i < 256; i++){for (int j = i; j < 256; j++){table2[i][j] = abs(i - j);table2[j][i] = table2[i][j];}}
}float MSE_block(Mat A, Mat B)
{float sum = 0.0;for (int j = 0; j < A.rows; j++){uchar *data1 = A.ptr<uchar>(j);uchar *data2 = B.ptr<uchar>(j);for (int i = 0; i < A.cols; i++){sum += table1[table2[data1[i]][data2[i]]];}}sum = sum / (A.rows*B.cols);return sum;
}//sigma越大越平滑
//halfKernelSize邻域窗
//halfSearchSize搜索窗
void NL_mean(Mat src, Mat &dst, double sigma, int halfKernelSize, int halfSearchSize)
{Mat boardSrc;dst.create(src.rows, src.cols, CV_8UC1);int boardSize = halfKernelSize + halfSearchSize;copyMakeBorder(src, boardSrc, boardSize, boardSize, boardSize, boardSize, BORDER_REFLECT);   //边界扩展double sigma_square = sigma*sigma;int rows = src.rows;int cols = src.cols;cal_lookup_table1();cal_lookup_table2();for (int j = boardSize; j < boardSize + rows; j++){uchar *dst_p = dst.ptr<uchar>(j - boardSize);for (int i = boardSize; i < boardSize + cols; i++){Mat patchA = boardSrc(Range(j - halfKernelSize, j + halfKernelSize), Range(i - halfKernelSize, i + halfKernelSize));double w = 0;double p = 0;double sumw = 0;for (int sr = -halfSearchSize; sr <= halfSearchSize; sr++)   //在搜索框内滑动{uchar *boardSrc_p = boardSrc.ptr<uchar>(j + sr);for (int sc = -halfSearchSize; sc <= halfSearchSize; sc++){Mat patchB = boardSrc(Range(j + sr - halfKernelSize, j + sr + halfKernelSize), Range(i + sc - halfKernelSize, i + sc + halfKernelSize));float d2 = MSE_block(patchA, patchB);w = exp(-d2 / sigma_square);p += boardSrc_p[i + sc] * w;sumw += w;}}dst_p[i - boardSize] = saturate_cast<uchar>(p / sumw);}}
}

参考文献

[1] 非局部均值滤波算法(NL-means)
[2] https://developer.aliyun.com/article/1156618
[3] opencv源码修改与使用:fastNlMeansDenoisingMulti()对CV_16U的支持


http://www.ppmy.cn/news/1520272.html

相关文章

WBGAPI:世界银行数据的现代化Python接口

WBGAPI:世界银行数据的现代化Python接口 项目简介 WBGAPI是一个为世界银行数据API提供现代、符合Python风格的访问工具。它旨在服务于数据新手和数据科学家。 特点 WBGAPI在几个方面与其他世界银行数据包不同,它充分利用了世界银行强大的API,同时减轻了一些特殊性的影响…

【个人笔记】python

Title: Python学习 1 基本用法 输入输出 格式化输出 python的格式化输出 str1 "num: {}".format(123) str1 "num: {:<44}".format(123) str1 "num: {:^20}".format(123){}与.format()&#xff1a;占位符、格式控制语句&#xff08;替换…

Java-树形图工具类TreeUtil

TreeUtil 工具类,包括列表转树形结构、遍历、查找和删除节点等功能。 import java.util.*;public class TreeUtil {/*** 将列表转换为树形结构。** @param target 扁平化的节点列表* @param getId 获取节点ID的函数* @param getParentId 获取节点父ID的函数* @…

kotlin中常用扩展函数

一、apply apply 用于在某个对象上执行一系列操作&#xff0c;然后返回该对象本身。这通常用于链式调用&#xff0c;以简化代码。 它被定义在 Any 类上&#xff0c;因此可以被任何类型的对象所使用。它接收一个 lambda 表达式作为参数&#xff0c;在这个 lambda 表达式中&#…

tcp/udp 可视化 调试工具; tcp/udp 发送客户端;查看tcp连接;netassist;packet sender;tcp view;

前言 最近对接物联设备&#xff0c;物联设备会往java端传送数据&#xff0c;在联调前&#xff0c;问厂商要来了样例数据&#xff0c;然后就找到了netassist和packetsender两款tcp发送工具&#xff0c;来模拟请求&#xff1b;还用tcpview来查看建立的连接。 netassist netass…

uniapp css样式盒子上面的箭头样式

在uniapp中&#xff0c;要创建一个盒子上面的箭头样式&#xff0c;可以使用CSS的伪元素::after或::before来实现。以下是一个简单的例子&#xff0c;展示了如何使用CSS创建一个向上的箭头。 <template><view class"box"><view class"box-conten…

一起学习LeetCode热题100道(60/100)

60.单词搜索(学习) 给定一个 m x n 二维字符网格 board 和一个字符串单词 word 。如果 word 存在于网格中&#xff0c;返回 true &#xff1b;否则&#xff0c;返回 false 。 单词必须按照字母顺序&#xff0c;通过相邻的单元格内的字母构成&#xff0c;其中“相邻”单元格是那…

Java-数据结构-链表-LinkedList(二)|ू・ω・` )

文本目录&#xff1a; ☛一、LinkedList&#xff08;无头双向非循环链表的结构&#xff09;&#xff1a; ☛ 二、LinkedList的模拟实现&#xff1a; ❄️1、操作方法&#xff1a; ▶&#xff08;1&#xff09;、display()方法&#xff1a; ▶&#xff08;2&#xff09;、size…