椒盐噪声是由图像传感器,传输信道,解码处理等产生的黑白相间的亮暗点噪声。椒盐噪声是指两种噪声,一种是盐噪声(salt noise)盐=白色(255),另一种是胡椒噪声(pepper noise),椒=黑色(0)。前者是高灰度噪声,后者属于低灰度噪声。一般两种噪声同时出现,呈现在图像上就是黑白杂点。对于彩色图像,也有可能表现为在单个像素BGR三个通道随机出现的255或0。椒盐噪声往往由图像切割引起,去除脉冲干扰及椒盐噪声最常用的算法是中值滤波。
Opencv提供了randu()函数生成单个均匀分布的随机数或随机数组。其函数原型为:
void randu(InputOutputArray dst,InputArray low,InputArray high); //参数说明: dst: //输出随机数组; 必须预先分配内存 low: //包含生成的随机数的下边界 high: //包含生成的随机数的上边界使用randu()函数给图像添加椒盐噪声的代码如下,pn是噪声密度(即包括噪声值的图像区域的百分比)。因此,大约有 p n × n u m e l ( f ) pn\times numel(f) pn×numel(f)个像素受到影响。默认的噪声密度为0.05。
void addSaltAndPepperNoise(InputArray srcImage, OutputArray dstImage, double pn=0.05) { Mat srcMat = srcImage.getMat(); srcImage.copyTo(dstImage); Mat dstMat = dstImage.getMat(); //图像通道判定 if (srcMat.channels() == 1) { Mat_<double> randux(srcMat.rows, srcMat.cols, CV_64F); randu(randux, Scalar_<double>(0.0), Scalar_<double>(1.0)); for(int i = 0; i < srcMat.rows; i++) { for(int j = 0; j < srcMat.cols; j++) { if(randux.at<double>(i, j) < pn / 2) { dstMat.at<uchar>(i, j) = 0; } else if((randux.at<double>(i, j) > pn / 2) & (randux.at<double>(i, j) < pn)) { dstMat.at<uchar>(i, j) = 255; } } } } else if(srcMat.channels() == 3) { Mat randux(srcMat.rows, srcMat.cols, CV_64FC3); randu(randux, Scalar(0.0, 0.0, 0.0), Scalar(1.0, 1.0, 1.0)); vector<Mat> dstchns; split(dstMat, dstchns); vector<Mat_<double>> xchns; split(randux, xchns); for(int k = 0; k < dstMat.channels(); k++) { for(int i = 0; i < dstMat.rows; i++) { for(int j = 0; j < dstMat.cols; j++) { if(xchns[k].at<double>(i, j) < pn / 2) { dstchns[k].at<uchar>(i, j) = 0; } else if((xchns[k].at<double>(i, j) > pn / 2) & (xchns[k].at<double>(i, j) < pn)) { dstchns[k].at<uchar>(i, j) = 255; } } } } merge(dstchns, dstMat); } } int main() { Mat srcImage = imread("lena.bmp", CV_LOAD_IMAGE_UNCHANGED); Mat noiseImage; addSaltAndPepperNoise(srcImage, noiseImage, 0.05); namedWindow("噪声图像"); imshow ("噪声图像", noiseImage); waitKey(0); return 0; }添加椒盐噪声后的图像为
所谓高斯噪声是指它的概率密度函数服从高斯分布(即正态分布)的一类加性噪声。 产生原因: 1)图像传感器在拍摄时不够明亮、亮度不够均匀; 2)电路各元器件自身噪声和相互影响; 3)图像传感器长期工作,温度过高。
我们可以使用C++ random库中的normal_distribution来产生正态分布随机数来添加高斯噪声。代码为:
void addGaussianNoise(InputArray _srcImage, OutputArray _dstImage, double mu = 0, double sigma = 0.1) { Mat srcImage = _srcImage.getMat(); _srcImage.copyTo (_dstImage); Mat _srcMat; int type = CV_MAKETYPE(CV_64F, srcImage.channels()); srcImage.convertTo (_srcMat, type, 1.0 / 255); Mat dstImage = _dstImage.getMat(); int channels = srcImage.channels(); int nRows = srcImage.rows; int nCols = srcImage.cols * channels; default_random_engine e; normal_distribution<double> n(mu, sigma); // 判断图像的连续性 if(srcImage.isContinuous()) { nCols *= nRows; nRows = 1; } for(int i = 0; i < nRows; ++i) { for(int j = 0; j < nCols; ++j) { // 添加高斯噪声 double val = _srcMat.ptr<double>(i)[j] + n(e); dstImage.ptr<uchar>(i)[j] = saturate_cast<uchar>(val * 255); } } }Box-Muller变换是通过服从均匀分布的随机变量,来构建服从高斯分布的随机变量的一种方法。基于这种变换,我们便可以得到一个从均匀分布中得到正态分布采样的算法。 定理:如果 U 1 U_1 U1和 U 2 U_2 U2是两个独立且服从(0,1)均匀分布的随机变量,则 X = − 2 ln U 1 cos ( 2 π U 2 ) X=\sqrt{-2 \ln U_1 }\cos(2\pi U_2) X=−2lnU1 cos(2πU2) Y = − 2 ln U 1 sin ( 2 π U 2 ) Y=\sqrt{-2 \ln U_1 }\sin(2\pi U_2) Y=−2lnU1 sin(2πU2) X X X, Y Y Y独立且服从标准正态分布(均值为0,方差为1)。 证明: 假设有两个标准的正态分布 X ∼ N ( 0 , 1 ) X \sim N(0,1) X∼N(0,1)和 Y ∼ N ( 0 , 1 ) Y\sim N(0,1) Y∼N(0,1),由于两者相互独立,则联合概率密度函数为: p ( x , y ) = p ( x ) p ( y ) = 1 2 π e − x 2 2 1 2 π e − x 2 2 = 1 2 π e − x 2 + y 2 2 p(x,y)=p(x)p(y)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2+y^2}{2}} p(x,y)=p(x)p(y)=2π 1e−2x22π 1e−2x2=2π 1e−2x2+y2 对上式进行极坐标变换,令 x = R cos θ x=R\cos\theta x=Rcosθ ,其中 θ ∼ [ 0 , 2 π ] θ\sim[0,2\pi] θ∼[0,2π],则有: F R ( r ) = F R ( R ≤ r ) = ∫ 0 2 π ∫ 0 r 1 2 π e − u 2 2 u d u d θ F_R(r)=F_R(R\leq r)=\int_0^{2\pi}\int_0^r\frac{1}{2\pi }e^{-\frac{u^2}{2}}udud\theta FR(r)=FR(R≤r)=∫02π∫0r2π1e−2u2ududθ 这个CDF函数的反函数可以写成 R = F − 1 ( u ) = − 2 ln ( 1 − u ) R=F^{-1}(u)=\sqrt{-2\ln(1-u)} R=F−1(u)=−2ln(1−u) 根据逆变换采样的原理,如果我们有个PDF为 p ( R ) p(R) p(R)的分布,那么对其CDF的反函数进行均匀采样所得的样本分布将符合 p ( R ) p(R) p(R)的分布,而如果u是均匀分布的,那么 U 1 = 1 − u U_1= 1-u U1=1−u也将是均匀分布的。于是用 U 1 U_1 U1 替换 1 − u 1-u 1−u,令 θ = 2 π U 2 \theta=2πU_2 θ=2πU2,其中 U 1 ∼ U ( 0 , 1 ) , U 2 ∼ U ( 0 , 1 ) U1\sim~U(0,1),U2\sim U(0,1) U1∼ U(0,1),U2∼U(0,1)最后可得: Z 0 = R cos θ = − 2 ln U 1 cos ( 2 π U 2 ) Z_0=R\cos\theta=\sqrt{-2\ln U_1}\cos(2\pi U_2) Z0=Rcosθ=−2lnU1 cos(2πU2) Z 1 = R cos θ = − 2 ln U 1 sin ( 2 π U 2 ) Z_1=R\cos\theta=\sqrt{-2\ln U_1}\sin(2\pi U_2) Z1=Rcosθ=−2lnU1 sin(2πU2) 因此,只需要有两个服从均匀分布的随机变量U1,U2,就可以得到一个服从正态分布的随机变量。 由概率论可知,当随机变量 X X X服从高斯分布,记为 X ∼ N ( μ 1 , σ 1 2 ) X\sim N(μ_1,σ_1^2) X∼N(μ1,σ12),那么 Y = A X + B Y=AX+B Y=AX+B也服从高斯分布,记为 Y ∼ N ( μ 2 , σ 2 2 ) Y\sim N(μ_2,σ_2^2) Y∼N(μ2,σ22),那么: μ 2 = E ( Y ) = E ( A X + B ) = A μ 1 + B μ_2=E(Y)=E(AX+B)=Aμ_1+B μ2=E(Y)=E(AX+B)=Aμ1+B σ 2 2 = D ( Y ) = D ( A X + B ) = A 2 D ( X ) = A 2 σ 1 2 σ_2^2=D(Y)=D(AX+B)=A^2D(X)=A^2σ_1^2 σ22=D(Y)=D(AX+B)=A2D(X)=A2σ12 当 μ 1 = 0 μ_1=0 μ1=0, σ 1 = 1 σ_1=1 σ1=1时, A = σ 2 , B = μ 2 A=σ_2,B=μ_2 A=σ2,B=μ2。由此可得,已知标准高斯分布 X ∼ N ( 0 , 1 ) X\sim N(0,1) X∼N(0,1),想要得到 Y ∼ N ( μ , σ 2 ) Y\sim N(μ,σ^2) Y∼N(μ,σ2)时,只需令 Y = σ X + μ Y=σX+μ Y=σX+μ
double generateGaussianNoise(double mu, double sigma) { // 定义小值 const double epsilon = std::numeric_limits<double>::min(); static double z0, z1; static bool flag = false; flag = !flag; // flag为假构造高斯随机变量X if (!flag) return z1 * sigma + mu; double u1, u2; // 构造随机变量 do { u1 = rand() * (1.0 / RAND_MAX); u2 = rand() * (1.0 / RAND_MAX); } while( u1 <= epsilon ); // flag为真构造高斯随机变量X z0 = sqrt(-2.0 * log(u1)) * cos(2 * CV_PI * u2); z1 = sqrt(-2.0 * log(u1)) * sin(2 * CV_PI * u2); return z0 * sigma + mu; } // 图像添加高斯噪声 void addGaussianNoise(InputArray _srcImage, OutputArray _dstImage, double mu = 0, double sigma = 0.1) { Mat srcImage = _srcImage.getMat(); _srcImage.copyTo (_dstImage); Mat _srcMat; int type = CV_MAKETYPE(CV_64F, srcImage.channels()); srcImage.convertTo (_srcMat, type, 1.0 / 255); Mat dstImage = _dstImage.getMat(); int channels = srcImage.channels(); int nRows = srcImage.rows; int nCols = srcImage.cols * channels; // 判断图像的连续性 if(srcImage.isContinuous()) { nCols *= nRows; nRows = 1; } for(int i = 0; i < nRows; ++i) { for(int j = 0; j < nCols; ++j) { // 添加高斯噪声 double val = _srcMat.ptr<double>(i)[j] + generateGaussianNoise(mu, sigma); dstImage.ptr<uchar>(i)[j] = saturate_cast<uchar>(val * 255); } } } int main() { Mat srcImage = imread("lena512.bmp"); Mat noiseImage; addGaussianNoise (srcImage, noiseImage); namedWindow("噪声图像"); imshow ("噪声图像", noiseImage); waitKey(0); return 0; }添加噪声后的图像如下图所示:
Opencv中的randn函数可以得到满足高斯分布的随机数,并填入数组或矩阵。函数原型为:
void randn (InputOutputArray dst, InputArray mean, InputArray stddev); dst;//输出数组或矩阵,必须预先分配内存,可以有1到4个通道 mean;//生成的随机数的平均值(期望值) stddev;//生成的随机数的标准差; 它可以是矢量(在这种情况下假设对角标准偏差矩阵)或方阵使用rand()函数得到添加高斯噪声的代码为:
void addGaussianNoise(InputArray srcImage, OutputArray dstImage, double mu = 0, double var = 0.1) { Mat srcMat = srcImage.getMat(); Mat _srcMat; int type = CV_MAKETYPE(CV_64F, srcMat.channels()); srcMat.convertTo (_srcMat, type, 1.0 / 255); Mat dstMat; if(srcMat.channels() == 1) { Mat randnx(srcMat.rows, srcMat.cols, CV_64F); randn(randnx, mu, var); add(_srcMat, randnx, dstMat, Mat(), CV_64F); } else if(srcMat.channels() == 3) { Mat randnx(srcMat.rows, srcMat.cols, CV_64FC3); randn(randnx, Scalar(mu, mu, mu), Scalar(var, var, var)); add(_srcMat, randnx, dstMat, Mat(), CV_64FC3); } int _type = CV_MAKETYPE(CV_8U, srcMat.channels()); dstMat.convertTo(dstMat, _type, 255); dstMat.copyTo(dstImage); } int main() { Mat srcImage = imread("lena.bmp"); Mat noiseImage; addGaussianNoise (srcImage, noiseImage, 0, 0.1); namedWindow("噪声图像"); imshow ("噪声图像", noiseImage); waitKey(0); return 0; }泊松噪声,就是符合泊松分布的噪声模型。由于光的量子特性,到达光电检测器表面的量子数目存在统计涨落,因此在低光度下检测图像存在颗粒性,这种颗粒性造成了图像对比度的变小和细节的遮掩。在低光度下,光子的发射可用泊松过程来描述,由它造成的测量不确定度称为泊松噪声。由光电检测器得到的图像信号包含在阵列单元所记录的“光电事件”平均数目中,这些“事件”的泊松分布意味着泊松噪声是一种与信号有关的噪声。泊松噪声既不是加性噪声,也不是乘性噪声,而是一种信号依赖噪声。对于一张图像而言,每个像素点的值都满足泊松分布,且每个像素点的泊松分布的均值是无噪图像在该像素点对应的值。 泊松分布噪声的概率质量函数为: P ( X = k ) = e − λ λ k k ! P(X=k)=\frac{e^{-\lambda}λ^k}{k!} P(X=k)=k!e−λλk 参数是λ单位时间(单位面积)内随机事件的平均发生率。服从泊松分布的随机变量,其数学期望与方差相等,同为参数λ: E ( X ) = D ( X ) = λ E(X)=D(X)=λ E(X)=D(X)=λ
一个用来生成随机泊松分布的数字(伪随机数抽样)的简单算法,已经由高德纳给出(见维基百科:泊松分布): 对于较大的λ值, e − λ e^{-λ} e−λ可能导致数值稳定性问题。对于较大λ值的一种解决方案是拒绝采样,另一种是采用泊松分布的高斯近似。
void addPoissonNoise(InputArray src, OutputArray dst) { Mat srcMat = src.getMat(); src.copyTo (dst); Mat dstMat = dst.getMat(); int channels = srcMat.channels(); int nRows = srcMat.rows; int nCols = srcMat.cols * channels; // 判断图像的连续性 if(srcMat.isContinuous()) { nCols *= nRows; nRows = 1; } default_random_engine e; uniform_real_distribution<double> u(0, 1); normal_distribution<double> n(0, 1); for(int i = 0; i < nRows; ++i) { for(int j = 0; j < nCols; ++j) { if(srcMat.ptr<uchar>(i)[j] <= 50) { long double L = exp(-srcMat.ptr<uchar>(i)[j]); int k = 0; long double p = 1.0; while(p > L) { p = p * u(e); k++; } dstMat.ptr<uchar>(i)[j] = k - 1; } else { //λ>50,采用高斯分布近似 dstMat.ptr<uchar>(i)[j] = saturate_cast<uchar>(srcMat.ptr<uchar>(i)[j] + sqrt(srcMat.ptr<uchar>(i)[j]) * n(e)); } } } }另外一个较为简单的方式是采用C++random库中的poisson_distribution生成服从泊松分布的随机数。代码如下:
void addPoissonNoise(InputArray src, OutputArray dst) { Mat srcMat = src.getMat(); src.copyTo (dst); Mat dstMat = dst.getMat(); int channels = srcMat.channels(); int nRows = srcMat.rows; int nCols = srcMat.cols * channels; // 判断图像的连续性 if(srcMat.isContinuous()) { nCols *= nRows; nRows = 1; } default_random_engine e; for(int i = 0; i < nRows; ++i) { for(int j = 0; j < nCols; ++j) { poisson_distribution<int> p(srcMat.ptr<uchar>(i)[j]); dstMat.ptr<uchar>(i)[j] = saturate_cast<uchar>(p(e)); } } } int main() { Mat srcImage = imread("lena512.bmp", CV_LOAD_IMAGE_UNCHANGED); Mat noiseImage; addPoissonNoise(srcImage, noiseImage); namedWindow("噪声图像"); imshow ("噪声图像", noiseImage); waitKey(0); return 0; }加入泊松噪声以后的图像为:
为了分析处理方便,往往将乘性噪声近似认为是加性噪声,而且总是假定信号和噪声是互相独立的。 为图像添加乘性噪声的代码如下,这里借鉴了Matlab中的imnoise中的处理方法。将乘性噪声添加到图像src上,其中噪声是均值为0,方差为var的均匀分布的随机噪声,var的默认值是0.04。
void addSpeckleNoise(InputArray src, OutputArray dst, double var = 0.04) { Mat srcMat = src.getMat(); Mat srcImage; Mat dstMat = dst.getMat(); int type = CV_MAKETYPE(CV_64F, srcMat.channels()); srcMat.convertTo (srcImage, type, 1.0 / 255); Mat dstImage = srcImage.clone(); if (srcMat.channels() == 1) { Mat_<double> randux(srcMat.rows, srcMat.cols, CV_64F); randu(randux, Scalar_<double>(0.0), Scalar_<double>(1.0)); randux = randux - 0.5; dstImage = srcImage * sqrt(12 * var); dstImage = srcImage + dstImage.mul(randux); } else if(srcMat.channels() == 3) { Mat randux(srcMat.rows, srcMat.cols, CV_64FC3); randu(randux, Scalar(0.0, 0.0, 0.0), Scalar(1.0, 1.0, 1.0)); randux = randux - Scalar(0.5, 0.5, 0.5); dstImage = srcImage * sqrt(12 * var); dstImage = srcImage + dstImage.mul(randux); } dstImage.convertTo(dst, srcMat.type(), 255); } int main() { Mat srcImage = imread("lena.bmp", CV_LOAD_IMAGE_UNCHANGED); Mat noiseImage; addSpeckleNoise(srcImage, noiseImage); namedWindow("噪声图像"); imshow ("噪声图像", noiseImage); waitKey(0); return 0; }