opencv-9-影像噪聲以及評估指標 PSNR 與SSIM

SChen1024發表於2020-05-01

開始之前

我們在將 opencv 的影像顯示在了 qt 的label 上, 我們能夠將圖顯示在label 上, 用於顯示我們的演算法,
我們在 opencv 上一篇文章中介紹了 opencv 的核操作, 我們這裡就要進入一個很重要的章節了,影像濾波操作, 也是影像核操作應用的一個很重要的章節,

那我們就從降噪的角度完整的講一下, 並通過 opencv 核的方式進行影像演算法操作, 【技術綜述】一文道盡傳統影像降噪方法 這篇文章寫的還算比較完整, 也是傳統的演算法的一個綜述過程,

目錄

前言

數字成像過程中由於電噪聲以及其他因素, 導致我們獲取到的影像存在噪聲,噪聲出現在輸入部分, 在後續的每個步驟都會受到影響, 所以在數字影像處理的前面必須要進行的一個步驟就是 影像降噪

每個做訊號處理的都會接觸到一類問題 , 訊號降噪, 讓人最頭疼的一門課,真是感謝老師給過, 但是後面自己用到的時候反而感覺真的好用, 原來是這樣, 然後就慢慢學會了怎麼使用吧..(感覺還是弱雞)
知乎可復現的影像降噪演算法總結這篇文章列出了一個能夠復現的影像降噪操作演算法列表, 近年來實現了的演算法可以見reproducible-image-denoising-state-of-the-art, 之後使用相應的文章進行演算法實現吧.( 又立了一個 flag )

影像降噪理論基礎

影像降噪主要的目的是在進行去除影像噪聲的同時保留儘可能多的主要特徵, 對於人眼來說, 區分噪聲還算比較容易, 但是對於計算機來說,輸入的都是資料, 我怎麼區分哪個是噪聲, 哪個不是噪聲呢, 這裡就要引入噪聲的理論基礎了

影像噪聲的產生

我們在之前的章節介紹了影像的程式系統, 實際上在成像過程中可能由於點噪聲, 量化過程等造成噪聲,
實際上的噪聲主要分為三種:

  1. 加性噪聲: 與輸入無關, \(f(x,y) = g(x,y) +n(x,y)\)
  2. 乘性噪聲: 與輸入訊號有關, \(f(x,y) = g(x,y) + n(x,y) \cdot g(x,y)\)
  3. 量化噪聲: 與輸入無關, 在影像量化過程噪聲的量化誤差導致的噪聲,

實際上後兩種很難解決, 目前處理的都是以 加性噪聲為主, 屬於隨機的噪聲訊號, 根據統計學的觀點, 噪聲在無限長時間窗內的噪聲和為0, 在第一類中的 \(n(x,y)\) 隨著時間存在正負訊號的不確定變化.

噪聲零和特點

上圖所示虛線代表真實訊號,紅藍線表示的就是隨機噪聲訊號,所有的隨機噪聲訊號求和後結果為0。

這裡關於噪聲的說明可以參考影像去噪演算法簡介

噪聲在理論上可以定義為“不可預測,只能用概率統計方法來認識的隨機誤差”。因此將影像噪聲看成是多維隨機過程是合適的,因而描述噪聲的方法完全可以借用隨機過程的描述,即用其概率分佈函式和概率密度分佈函式。但在很多情況下,這樣的描述方法是很複雜的,甚至是不可能的。而實際應用往往也不必要。通常是用其數字特徵,即均值方差,相關函式等。因為這些數字特徵都可以從某些方面反映出噪聲的特徵。

我認為影像噪聲的成因分類與常見影像去噪演算法簡介這篇文章關於噪聲的分類部分講的還比較細, 可以參考

影像噪聲的模型

由於我們認為噪聲在時間尺度的隨機性, 但是我們可以使用噪聲的概率分佈與概率密度函式進行描述, 那麼我們就能將噪聲根據其分佈特點進行分類,
我們稍微介紹一下常見的噪聲模型吧

一些重要的噪聲概率密度函式

噪聲模型主要可以分為:

  • 高斯噪聲,高斯噪聲模型經常被用於實踐中。
  • 脈衝噪聲(椒鹽噪聲),影像上一個個點,也可稱為散粒和尖峰噪聲。
  • 伽馬噪聲
  • 瑞利噪聲
  • 指數分佈噪聲
  • 均勻分佈噪聲

這裡能查到的資料很多, 可以看我們的參考部分, 內容都一樣, 再寫只是浪費時間和精力, 有興趣的可以自己翻閱

影像降噪操作

其實吧, 我就不應該講那麼多, 直接開始影像處理部分就行了, 為了開始進行影像處理, 我們要先進行一點小工作, 我們要按造以下步驟進行降噪演算法的比較,

  1. 選擇標準影像--- lena.png
  2. 新增噪聲
  3. 量化噪聲
  4. 降噪操作
  5. 量化結果值
  6. 比較結果

在我們進行演算法比對之前, 我們選擇的是 lena 的影像, 加入隨機噪聲, 然後計算出來 一個噪聲的比例, 進行降噪操作, 再次計算以下噪聲引數, 看下效果值.

如果是進行演算法比較的時候, 最好選擇現有的降噪的資料集進行比較, 比如, Kodak , BSD

噪聲新增

我們認為噪聲是隨機的, 我們生成隨機數加在原始影像上便能夠得到噪聲影像, opencv 沒有提供相應的實現, 但是知道原理了, 寫起來都比較簡單, 我比較喜歡
影像處理基礎(1):噪聲的新增和過濾 使用的方法, 他使用的是 梅森旋轉演算法 來實現的偽隨機演算法,

其實吧這裡我也不懂, 但是隨機數能用就行了, 我又不是數學家, 然後看到了 談談梅森旋轉:演算法及其爆破

這裡就不重複造輪子了, 直接複製他給出的程式碼就好,

// 新增椒鹽噪聲 // 生成 隨機 num 個 白點
void addSaltNoise(Mat &m, int num)
{
	// 隨機數產生器
	std::random_device rd; //種子
	std::mt19937 gen(rd()); // 隨機數引擎

	auto cols = m.cols * m.channels();

	for (int i = 0; i < num; i++)
	{
		auto row = static_cast<int>(gen() % m.rows);
		auto col = static_cast<int>(gen() % cols);

		auto p = m.ptr<uchar>(row);
		p[col++] = 255;
		p[col++] = 255;
		p[col] = 255;
	}
}

// 新增Gussia噪聲
// 使用指標訪問
void addGaussianNoise(Mat &m, int mu, int sigma)
{
	// 產生高斯分佈隨機數發生器
	std::random_device rd;
	std::mt19937 gen(rd());

	std::normal_distribution<> d(mu, sigma);

	auto rows = m.rows; // 行數
	auto cols = m.cols * m.channels(); // 列數

	for (int i = 0; i < rows; i++)
	{
		auto p = m.ptr<uchar>(i); // 取得行首指標
		for (int j = 0; j < cols; j++)
		{
			auto tmp = p[j] + d(gen);
			tmp = tmp > 255 ? 255 : tmp;
			tmp = tmp < 0 ? 0 : tmp;
			p[j] = tmp;
		}
	}
}

噪聲量化方法

這裡其實涉及到影像質量評估的領域,可以參考影像質量評價概述(評估指標、傳統檢測方法)介紹的方法, 存在太多的計算方式,

我們必須選擇一個量化噪聲的方式進行影像質量的評估, 一般進行噪聲評估手段就是噪聲比(Signal to Noise Ratio,SNR),峰值訊雜比(Peak Signal to Noise Ratio, PSNR) , 均方差值(Mean Square Error, MSE), 結構相似性(Structural SIMilarity, SSIM),

我們一個一個來看, 均方差值是用於比較兩幅影像 K, I 的均方差值

\[MSE=\frac{1}{mn}\sum_{i=0}^{n-1}\sum_{j=0}^{m-1}\parallel K(i,j)-I(i,j)\parallel^{2} \]

峰值訊雜比PSNR衡量影像失真或是噪聲水平的客觀標準。2個影像之間PSNR值越大,則越相似。普遍基準為30dB,30dB以下的影像劣化較為明顯。定義為,

\[PSNR=10log_{10}(\frac{MAX^{2}}{MSE}) \]

其中\(MAX^2\) 為圖片可能的最大畫素值。如果每個畫素都由 8 位二進位制來表示,那麼就為 255。

SNR用於描述訊號與噪聲的比值

\[SNR (dB)=10 log_{10} \left[ \frac{\sum_{x=0}^{m-1} \sum_{y=0}^{n-1}(f(x, y))^{2}}{\sum_{x=0}^{m-1} \sum_{y=0}^{n-1}(f(x, y)-\hat{f}(x, y))^{2}}\right] \]

SSIM 描述兩個影像的相似性, 通過三個進行比較, 亮度,對比度和結構, 參考影像質量評價指標之 PSNR 和 SSIM

\[l(x, y)=\frac{2 \mu_{x} \mu_{y}+c_{1}}{\mu_{x}^{2}+\mu_{y}^{2}+c_{1}} c(x, y)=\frac{2 \sigma_{x} \sigma_{y}+c_{2}}{\sigma_{x}^{2}+\sigma_{y}^{2}+c_{2}} s(x, y)=\frac{\sigma_{x y}+c_{3}}{\sigma_{x} \sigma_{y}+c_{3}} \]

\[\operatorname{SSIM}(x, y)=\frac{\left(2 \mu_{x} \mu_{y}+c_{1}\right)\left(2 \sigma_{x y}+c_{2}\right)}{\left(\mu_{x}^{2}+\mu_{y}^{2}+c_{1}\right)\left(\sigma_{x}^{2}+\sigma_{y}^{2}+c_{2}\right)} \]

一般取\(c_3 = \frac{c_2}{2}\)
\(u_x\)\(x\) 的均值
\(u_y\)\(y\) 的均值
\(\sigma_x^2\)\(x\) 的方差
\(\sigma_y^2\)\(y\) 的方差
\(\sigma_{xy}\)\(x\)\(y\) 的協方差
\(c_1 = (k_1 L)^2, c_2=(k_2 L)^2\) 為兩個常數,避免除零
\(L\) 為畫素值的範圍,\((0,255)\)
\(k_1 = 0.01, k_2 = 0.03\) 為預設值
預設引數\(\alpha = 1, \beta = 1, \gamma = 1\)

opencv 計算 PSNR 和 SSIM

本來不想寫這麼多的, 但是 opencv 給出了一個例程Similarity check (PNSR and SSIM) on the GPU, 提供了計算的方法, 自己不用去寫了, 豈不是很爽, 所以上面就詳細介紹了各個方法的使用.
官方給出了普通版本以及 GPU 加速的版本, 我們暫時只使用基礎的版本就好,
PSNR返回一個浮點數,如果兩個輸入在30到50之間相似(越高越好)。
SSIM返回影像的MSSIM。這也是一個介於零和一之間的浮點數(越高越好),但是每個通道都有一個浮點數。因此,我們返回一個Scalar OpenCV資料結構:

double getPSNR(const Mat& I1, const Mat& I2)
{
    Mat s1;
    absdiff(I1, I2, s1);       // |I1 - I2|
    s1.convertTo(s1, CV_32F);  // cannot make a square on 8 bits
    s1 = s1.mul(s1);           // |I1 - I2|^2
    Scalar s = sum(s1);         // sum elements per channel
    double sse = s.val[0] + s.val[1] + s.val[2]; // sum channels
    if( sse <= 1e-10) // for small values return zero
        return 0;
    else
    {
        double  mse =sse /(double)(I1.channels() * I1.total());
        double psnr = 10.0*log10((255*255)/mse);
        return psnr;
    }
}
Scalar getMSSIM( const Mat& i1, const Mat& i2)
{
    const double C1 = 6.5025, C2 = 58.5225;
    /***************************** INITS **********************************/
    int d     = CV_32F;
    Mat I1, I2;
    i1.convertTo(I1, d);           // cannot calculate on one byte large values
    i2.convertTo(I2, d);
    Mat I2_2   = I2.mul(I2);        // I2^2
    Mat I1_2   = I1.mul(I1);        // I1^2
    Mat I1_I2  = I1.mul(I2);        // I1 * I2
    /*************************** END INITS **********************************/
    Mat mu1, mu2;   // PRELIMINARY COMPUTING
    GaussianBlur(I1, mu1, Size(11, 11), 1.5);
    GaussianBlur(I2, mu2, Size(11, 11), 1.5);
    Mat mu1_2   =   mu1.mul(mu1);
    Mat mu2_2   =   mu2.mul(mu2);
    Mat mu1_mu2 =   mu1.mul(mu2);
    Mat sigma1_2, sigma2_2, sigma12;
    GaussianBlur(I1_2, sigma1_2, Size(11, 11), 1.5);
    sigma1_2 -= mu1_2;
    GaussianBlur(I2_2, sigma2_2, Size(11, 11), 1.5);
    sigma2_2 -= mu2_2;
    GaussianBlur(I1_I2, sigma12, Size(11, 11), 1.5);
    sigma12 -= mu1_mu2;
    Mat t1, t2, t3;
    t1 = 2 * mu1_mu2 + C1;
    t2 = 2 * sigma12 + C2;
    t3 = t1.mul(t2);              // t3 = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))
    t1 = mu1_2 + mu2_2 + C1;
    t2 = sigma1_2 + sigma2_2 + C2;
    t1 = t1.mul(t2);               // t1 =((mu1_2 + mu2_2 + C1).*(sigma1_2 + sigma2_2 + C2))
    Mat ssim_map;
    divide(t3, t1, ssim_map);      // ssim_map =  t3./t1;
    Scalar mssim = mean( ssim_map ); // mssim = average of ssim map
    return mssim;
}

演算法噪聲資料

我們完成了噪聲新增以及噪聲的量化, 我們來試一下, 給影像隨機新增一定的噪聲, 然後看下相應的引數變化情況對比來看就好

椒鹽噪聲測試

我們先來測試椒鹽噪聲 分別計算沒有噪聲的圖, 以及新增了 1000個 和10000個噪聲的資料結果, 並將後面兩個顯示出來


void MainWindow::testFunc1(void)
{
    // 新增椒鹽噪聲 並計算 PSNR和 SSIM
    cv::Mat salt_img;
    double psnr = 0;
    cv::Scalar mssim;

    QString res_temp = "Salt-%1 : psnr:%2, mssim: B:%3 G:%4 R:%5 ";
    QString res_str;

    // 計算三組影像的引數 0, 1000, 10000

    // 複製原始影像, 新增噪聲, 計算 psnr和ssim  顯示在 ui上
    salt_img = gSrcImg.clone();
    addSaltNoise(salt_img,0);

    psnr = getPSNR(gSrcImg, salt_img);
    mssim = getMSSIM(gSrcImg,salt_img);
    res_str = res_temp.arg(0)
                        .arg(psnr)
                        .arg(mssim.val[0])
                        .arg(mssim.val[1])
                        .arg(mssim.val[2]);
    ui->pt_log->appendPlainText(res_str);

    salt_img = gSrcImg.clone();
    addSaltNoise(salt_img,1000);

    psnr = getPSNR(gSrcImg, salt_img);
    mssim = getMSSIM(gSrcImg,salt_img);
    res_str = res_temp.arg(1000)
                        .arg(psnr)
                        .arg(mssim.val[0])
                        .arg(mssim.val[1])
                        .arg(mssim.val[2]);
    ui->pt_log->appendPlainText(res_str);

    // 左側顯示 1000 噪聲 右側顯示 10000 噪聲
    ShowMatOnQtLabel(salt_img,ui->lb_src);

    salt_img = gSrcImg.clone();
    addSaltNoise(salt_img,10000);

    psnr = getPSNR(gSrcImg, salt_img);
    mssim = getMSSIM(gSrcImg,salt_img);
    res_str = res_temp.arg(10000)
                        .arg(psnr)
                        .arg(mssim.val[0])
                        .arg(mssim.val[1])
                        .arg(mssim.val[2]);
    ui->pt_log->appendPlainText(res_str);

    ShowMatOnQtLabel(salt_img,ui->lb_dst);
}

1000 和 10000 的椒鹽噪聲對比

我們可以直接計算得到椒鹽噪聲 psnr 和 ssim 都是越大越好的, 可以明顯的看到影像質量退化

Salt-0 : psnr:0, mssim: B:1 G:1 R:1 
Salt-1000 : psnr:27.7528, mssim: B:0.865341 G:0.870555 R:0.914122 
Salt-10000 : psnr:17.8062, mssim: B:0.311999 G:0.327485 R:0.493874 

高斯噪聲測試

高斯噪聲我們測試了四組 分別使用引數(0,1) (0,10)(10,1)(10,10) 作為高斯引數, 最終得到後面的圖, 然後計算得到的結果, 我們做的結果比較簡單, 可以參考數字影像處理——新增高斯噪聲&椒鹽噪聲, 給出了很多的圖, 可以參考學


void MainWindow::testFunc2(void)
{
    // 新增高斯噪聲 並計算 PSNR和 SSIM
    cv::Mat guass_img;
    double psnr = 0;
    cv::Scalar mssim;

    QString res_temp = "gauss-%1- %2 : psnr:%3, mssim: B:%4 G:%5 R:%6 ";
    QString res_str;

    // 計算三組影像的引數 (0,1) (0,10), (10,1), (10,10)

    // 複製原始影像, 新增噪聲, 計算 psnr和ssim  顯示在 ui上
    guass_img = gSrcImg.clone();
    addGaussianNoise(guass_img,0,1);

    psnr = getPSNR(gSrcImg, guass_img);
    mssim = getMSSIM(gSrcImg,guass_img);
    res_str = res_temp.arg(0)
                        .arg(1)
                        .arg(psnr)
                        .arg(mssim.val[0])
                        .arg(mssim.val[1])
                        .arg(mssim.val[2]);
    ui->pt_log->appendPlainText(res_str);
    guass_img = gSrcImg.clone();
    addGaussianNoise(guass_img,0,10);

    psnr = getPSNR(gSrcImg, guass_img);
    mssim = getMSSIM(gSrcImg,guass_img);
    res_str = res_temp.arg(0)
                        .arg(10)
                        .arg(psnr)
                        .arg(mssim.val[0])
                        .arg(mssim.val[1])
                        .arg(mssim.val[2]);
    ui->pt_log->appendPlainText(res_str);

    guass_img = gSrcImg.clone();
    addGaussianNoise(guass_img,10,1);

    psnr = getPSNR(gSrcImg, guass_img);
    mssim = getMSSIM(gSrcImg,guass_img);
    res_str = res_temp.arg(10)
                        .arg(1)
                        .arg(psnr)
                        .arg(mssim.val[0])
                        .arg(mssim.val[1])
                        .arg(mssim.val[2]);
    ui->pt_log->appendPlainText(res_str);

    guass_img = gSrcImg.clone();
    addGaussianNoise(guass_img,10,10);

    psnr = getPSNR(gSrcImg, guass_img);
    mssim = getMSSIM(gSrcImg,guass_img);
    res_str = res_temp.arg(10)
                        .arg(10)
                        .arg(psnr)
                        .arg(mssim.val[0])
                        .arg(mssim.val[1])
                        .arg(mssim.val[2]);
    ui->pt_log->appendPlainText(res_str);
}

四組高斯測試結果對比

gauss-0- 1 : psnr:46.8791, mssim: B:0.991811 G:0.991622 R:0.992751 
gauss-0- 10 : psnr:28.1229, mssim: B:0.614219 G:0.608773 R:0.648285 
gauss-10- 1 : psnr:28.5293, mssim: B:0.978448 G:0.980308 R:0.987926 
gauss-10- 10 : psnr:25.3511, mssim: B:0.605665 G:0.600491 R:0.646768 

總結

原本想把濾波一起做了的, 但是越寫越, 就不做太多的處理了, 我們算是介紹了噪聲的來源, 噪聲的模型, 以及個噪聲的量化方式,
然後介紹了影像新增噪聲的方法 我們分別給影像新增椒鹽噪聲與高斯噪聲, 然後分別量化了噪聲的結果值, 進行對比展示,

示例的圖不是很多, 程式是在程式碼庫裡面的, 可以直接去自己實現, 然後進行 進行更多圖的展示

參考

  1. 《高斯噪聲_百度百科》. 見於 2020年4月30日. https://baike.baidu.com/item/高斯噪聲.
  2. 知乎專欄. 《【技術綜述】一文道盡傳統影像降噪方法》. 見於 2020年4月29日. https://zhuanlan.zhihu.com/p/51403693.
  3. 知乎專欄. 《可復現的影像降噪演算法總結》. 見於 2020年4月29日. https://zhuanlan.zhihu.com/p/32502816.
  4. 《梅森旋轉演算法》. 收入 維基百科,自由的百科全書, 2019年11月4日. https://zh.wikipedia.org/w/index.php?title=梅森旋轉演算法&oldid=56745942.
  5. 《實現灰度影像峰值訊雜比計算_人工智慧_松子茶的專欄-CSDN部落格》. 見於 2020年4月30日. https://blog.csdn.net/songzitea/article/details/17529445.
  6. 《數字影像處理-噪聲 - Mohanson》. 見於 2020年4月30日. http://accu.cc/content/pil/noise/.
  7. 《影像處理基礎(1):噪聲的新增和過濾 - Brook_icv - 部落格園》. 見於 2020年4月30日. https://www.cnblogs.com/wangguchangqing/p/6372025.html.
  8. 《影像處理PSNR及其計算(OpenCV和matlab實現)_人工智慧_無機器不學習-加大碼的分享-CSDN部落格》. 見於 2020年4月30日. https://blog.csdn.net/laoxuan2011/article/details/51519062.
  9. 《影像的 SNR 和 PSNR 的計算 - rldts - 部落格園》. 見於 2020年4月30日. https://www.cnblogs.com/qrlozte/p/5340216.html.
  10. 《影像去噪演算法簡介 - InfantSorrow - 部落格園》. 見於 2020年4月29日. https://www.cnblogs.com/CCBB/archive/2011/01/06/1929033.html.
  11. 《影像噪聲的成因分類與常見影像去噪演算法簡介_Java_qq_27606639的部落格-CSDN部落格》. 見於 2020年4月30日. https://blog.csdn.net/qq_27606639/article/details/80912071.
  12. 《影像質量評估指標 SSIM / PSNR / MSE_人工智慧_兔角與禪-CSDN部落格》. 見於 2020年4月30日. https://blog.csdn.net/edogawachia/article/details/78756680.
  13. 《影像質量評價概述(評估指標、傳統檢測方法)_人工智慧_qq_23304241的部落格-CSDN部落格》. 見於 2020年4月30日. https://blog.csdn.net/qq_23304241/article/details/80953613.
  14. 《影像降噪》. 收入 維基百科,自由的百科全書, 2018年9月20日. https://zh.wikipedia.org/w/index.php?title=影像降噪&oldid=51354600.

相關文章