技術標籤：矩陣的奇異值分解

注：本文是蔡樂衡同學對奇異值分解在影象去噪中應用的介紹

1 奇異值分解在影象去噪中的應用

將影象去噪問題用數學式一般化為：

其中訊號矩陣為X，噪聲為D，我們觀測到的含噪矩陣記為Y。我們的目標是從觀測到的含噪矩陣Y中恢復訊號矩陣X。

基於奇異值分解的去噪技術屬於子空間演算法的一種，簡單來說我們希望將含噪矩陣的向量空間分解為分別由純淨訊號主導的和噪聲訊號主導的兩個子空間，進而去除落在噪聲空間中的含噪訊號向量分量來估計純淨訊號。

類似奇異值分解在影象壓縮中的應用，一個非常普遍的辦法是：首先將影象對應的畫素矩陣進行奇異值分解，提取出其中較大的奇異值，我們認為這些奇異值代表了影象主要的訊號部分；截斷較小的奇異值，因為它們很可能代表了影象中某些噪聲，再將提取出的奇異值與對應特徵向量組合復原影象即可。

下圖為本文所用示例的不含噪音的訊號影象(大小為255*255的灰度影象)

2 稀疏的影象噪聲

本文中的稀疏噪聲表示在影象某些區域性存在的小汙點，我們通過在原始訊號矩陣X中隨機稀疏灑點，並在這些灑點位置處新增服從正態分佈的隨機擾動作為噪音D，即可得到包含稀疏噪聲的觀測影象Y：圖中的白色、黑色小點即為噪聲。

我們直接對觀測矩陣Y進行奇異值分解，截斷前25個奇異值後復原影象，即可得到去噪後的影象如下：儘管某些邊緣輪廓出現模糊，但幾乎所有的噪音都被有效消除了！

3 稠密的影象噪聲

當原始矩陣每一個畫素點都存在噪音時，奇異值分解仍然有效嗎？

3.1 基於低秩模型的方法

基於低秩模型方法和上述處理稀疏噪聲的方法相同，都是尋找F範數下秩為k的矩陣X的最優近似。即：我們直接對含噪矩陣Y進行奇異值分解，可以最小化如下平方誤差：

該方法求解的可表示為：

其中和分別為含噪矩陣Y的左右奇異向量，是Y的k個較大的奇異值。上圖展示了k取25時的去噪效果。

3.2 基於迴歸的方法

這裡給定含噪矩陣Y，我們需要尋找矩陣H以最小化如下誤差：

此時最優的H表示為：。但由於上式中X未知，僅通過Y我們不能得到。因此我們不妨做一些假設：

其中為噪聲方差(噪聲方差的估計請參考文獻【1】) 進而我們可以得到：

我們可以看到這兩種方法估計的具有相同的奇異向量，但具有不同的奇異值。上圖為方法二中k取25時的去噪影象。經比較，兩種方法提取的訊號影象幾乎相同，去噪效果遠不如對稀疏噪聲進行去噪。筆者認為一個可能的原因在於，當每一個畫素點都有噪音干擾時，那些較大的奇異值以及對應的奇異向量也受到了噪聲干擾，因此直接提取或直接復原影象的去噪效果並不理想。

目前學術上已經有較多的文獻在常規SVD方法基礎上進行改進，以更好地提取出訊號影象，如參考文獻2、3。

4 程式碼部分

library(wavelets)
library(jpeg)
library(magick)#載入magick包

setwd('C:/Users/micha/Pictures')
#圖片路徑
path='cat.jpg'
img#讀入圖片
#檢視影象屬性
result=image_info(img)
image_info(img)

#獲取圖片矩陣行數H列數W
H=result$height
W=result$width

#轉為灰度圖
gray'gray')
gray#顯示灰度影象

#轉為rgb圖
rgb'rgb')
rgb#顯示彩色影象

#獲取圖片灰度值矩陣X
X1]
#歸一化
X255
#儲存灰度圖
writeJPEG(X,'X.jpg')


#對原始影象奇異值分解
svd_result1=svd(X)
u1=svd_result1$u
v1=svd_result1$v
d1=svd_result1$d
d1=diag(d1)
#截斷奇異值k
k=25
svd_X=u1[,1:k]%*%d1[1:k,1:k]%*%t(v1[,1:k])
svd_X[svd_X<0]=0
svd_X[svd_X>1]=1
#儲存原始影象的壓縮影象
writeJPEG(svd_X,'svd_X.jpg')

#生成高斯噪聲
sigma=0.05
D0,sigma),H,W)
#合成含噪影象矩陣Y
YY[Y<0]=0
Y[Y>1]=1
#儲存含噪影象
writeJPEG(Y,'Y.jpg')

#對含噪矩陣直接奇異值分解
svd_result2=svd(Y)
u2=svd_result2$u
v2=svd_result2$v
d2=svd_result2$d
d2=diag(d2)
#截斷奇異值k
svd_Y=u2[,1:k]%*%d2[1:k,1:k]%*%t(v2[,1:k])
svd_Y[svd_Y<0]=0
svd_Y[svd_Y>1]=1
writeJPEG(svd_Y,'svd_Y.jpg')

#最小方差估計
u3=svd_result2$u
v3=svd_result2$v
d3=svd_result2$d

#多尺度小波分解
wt2)
wt_W=[email protected]$W1
#估計噪聲標準差
median(abs(wt_W-median(wt_W)))/qnorm(0.75,0,1)
sigmahat=median(abs(wt_W-median(wt_W)))/qnorm(0.75,0,1)

#對奇異值進行調整
d3=(d3^2-sigmahat^2)/d3
d3=diag(d3)
#截斷奇異值k
svd_Y1=u3[,1:k]%*%d3[1:k,1:k]%*%t(v3[,1:k])
svd_Y1[svd_Y1<0]=0
svd_Y1[svd_Y1>1]=1
writeJPEG(svd_Y1,'svd_Y1.jpg')

#############################
#稀疏噪音
#生成高斯噪聲
sigma=0.6
D0,H,W)
#噪聲點數目
num=200
temp1=sample(1:H,num,T)
temp2=sample(1:W,num,T)
for(iinseq(1,num)){
D[temp1[i],temp2[i]]=rnorm(1,0,0.3)
}

#合成含噪影象矩陣Y
YY[Y<0]=0
Y[Y>1]=1
#儲存含噪影象
writeJPEG(Y,'yy.jpg')

#對含噪矩陣直接奇異值分解
svd_result2=svd(Y)
u2=svd_result2$u
v2=svd_result2$v
d2=svd_result2$d
d2=diag(d2)
#截斷奇異值k
svd_Y=u2[,1:k]%*%d2[1:k,1:k]%*%t(v2[,1:k])
svd_Y[svd_Y<0]=0
svd_Y[svd_Y>1]=1
writeJPEG(svd_Y,'svd_yy.jpg')

5 參考文獻

【1】Donoho D L . De-noising by soft thresholding[J]. IEEE Transactions on Information Theory, 1995.

【2】Nadakuditi R R . OptShrink: An algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage[J]. IEEE Transactions on Information Theory, 2013, 60(5):3002-3018.

【3】Gavish, Matan, and Donoho, David L. "Optimal Shrinkage of Singular Values." IEEE Transactions on Information Theory 63.4 (2017): 2137-152. Web.