UA MATH575B 數值分析下III 圖像恢復

• 圖像去噪的優化模型
• 算法實現

圖像去噪的優化模型

假設$V$代表原圖，觀測到的圖像是$W$，$W$是原圖與噪聲的疊加$W=V+\xi$，我們的目標是對$W$去噪，恢復原來的圖像，可以使用一個優化模型來完成這個任務。假設恢復后的圖像為$U$，則
$$U^{?}=\underset{U}{\mathrm{argmin}}(\frac{1}{2}{\Vert U?W\Vert }_2^2+\lambda {\Vert ?U\Vert }_p^{\nu })$$
這個優化目標的第一項${\Vert U?W\Vert }_2^2$作用是保證恢復后的原圖與有噪聲的圖像差別不會太大；第二項的作用是通過系數$\lambda$來調節恢復后的圖像平滑程度，通常設定$\nu =p$。為了讓這個模型更可行，將圖像（image）想象成一個圖（graph），圖的頂點是圖像的像素點，圖的邊代表其兩個端點對應的像素點相鄰。從而優化目標的第一項可以表示為
$${\Vert U?W\Vert }_2^2=\sum _v(U_v?W_v{)}^2$$
$v$代表圖像的像素點。第二項可以表示為
$${\Vert ?U\Vert }_p^p=\sum _{e=v_1{v}_2}|U_{v_1}?U_{v_2}{|}^p$$
平滑項取$L_1$范數的效果一般是比較好的，所以令$p=1$，但這樣第二項就成了絕對值函數不好做優化，因此考慮對絕對值做一個近似
$$|U_{v_1}?U_{v_2}|\approx \sqrt{{?}^2+(U_{v_1}?U_{v_2}{)}^2}$$
其中$?$是一個絕對很小的常數。因此優化問題可以表示為
$$U^{?}=\underset{U}{\mathrm{argmin}}(\frac{1}{2}\sum _v(U_v?W_v{)}^2+\lambda \sum _{e=v_1{v}_2}\sqrt{{?}^2+(U_{v_1}?U_{v_2}{)}^2}）$$
現在優化目標是一個平滑的凸函數了，記為$f(U)$。其一階偏導為
$$\frac{?f}{?U_v}=U_v?W_v+\lambda \sum _{v^{\prime }}\frac{U_v?U_{v^{\prime }}}{\sqrt{{?}^2+(U_v?U_{v^{\prime }}{)}^2}}$$
$v^{\prime }$代表與$v$相鄰的像素點。對所有的像素點$v$求解一階偏導等于零的非線性方程可以得到恢復的圖像。但這里考慮采用梯度下降算法來做，但并不采用傳統的梯度下降每一步循環計算梯度搜索步長更新輸入的方式，這里嘗試用數值ODE方法來做梯度下降。令
$$\frac{d{U}_v}{dt}=?\frac{?f}{?U_v}=W_v?U_v?\lambda \sum _{v^{\prime }}\frac{U_v?U_{v^{\prime }}}{\sqrt{{?}^2+(U_v?U_{v^{\prime }}{)}^2}}$$
其中$t$是一個參數，使用數值算法求解這個ODE，當$U_v(t)$滿足停止準則后將其值輸出作為圖像在像素點$v$的取值。但這樣有一個壞處，如果梯度下降比較順利的話收斂時比較快的，但不考慮自適應步長的ODE算法（比如RK類算法）步長都是固定的，因此即使能比較順利的下降收斂也不會變快。為了改掉這個缺陷，可以仿照加速梯度下降（Nesterov’s accelerated gradient descent，AGD）方法，在ODE中引入一個慣性項$m\frac{d^2{U}_v}{d{t}^2}$，從而
$$m\frac{d^2{U}_v}{d{t}^2}+\frac{d{U}_v}{dt}=W_v?U_v?\lambda \sum _{v^{\prime }}\frac{U_v?U_{v^{\prime }}}{\sqrt{{?}^2+(U_v?U_{v^{\prime }}{)}^2}}$$
$m$的物理意義是質量，但在計算中它就是一個超參。

算法實現

為什么要用這個枯燥的例子呢？因為這是我的作業【吐血.jpg】下面這張灰度圖是我要用來去噪的，不難看出它居然是個A！

這張圖可以用一個矩陣來表示，對應上面的記號，用$W$來表示，大概長下面這樣，每個值代表那個像素點的灰度，取值是0-255，每個值占8 bits。上面的原圖就是用這個矩陣生成的，在matlab里可以用這個來生成

function outGrey(A) mat = A; %# A n-by-n matrix imagesc(mat); %# Create a colored plot of the matrix values colormap(flipud(gray)); end

這張圖是$64\times 64$的，這就意味著一共有4096個像素點，需要重構的圖像$U$是一個有4096個元素的矩陣。在這個ODE-based梯度下降的計算模型中
$$m\frac{d^2{U}_v}{d{t}^2}+\frac{d{U}_v}{dt}=W_v?U_v?\lambda \sum _{v^{\prime }}\frac{U_v?U_{v^{\prime }}}{\sqrt{{?}^2+(U_v?U_{v^{\prime }}{)}^2}}$$
參考老師給的參數，$m=0.2,?=0.1$。先把這個兩階的ODE用二元一階ODE表示一下，假設$X_v=\frac{d{U}_v}{dt}$，則這個方程可以寫成
$$\begin{gathered} \frac{d{X}_v}{dt}=(W_v?U_v?\lambda \sum _{v^{\prime }}\frac{U_v?U_{v^{\prime }}}{\sqrt{{?}^2+(U_v?U_{v^{\prime }}{)}^2}}?X_v)/m \\ \frac{d{U}_v}{dt}=X_v \end{gathered}$$
這個看起來是二元一階ODE的東西其實有8192個未知數。現在把這個圖像讀進Matlab，再定義一些參數

global W lambda epsilon mass; W = textread('A051.txt');% Read image from txt file epsilon = .1; mass = .2;% Set parameters

定義一下ODE右邊那個函數，RK4會用到，要注意四個角和四條邊上的像素相鄰像素點分別只有2和3個。

function f = RHS(U,X,i,j) %% This is the right hand side function of ODE % input i,j indicate the location of pixel global W lambda mass; S = 0; if ((i==1) &&(j==1))S = S + (U(i,j) - U(i+1,j))/sqrt(epsilon^2 + (U(i,j) - U(i+1,j))^2);S = S + (U(i,j) - U(i,j+1))/sqrt(epsilon^2 + (U(i,j) - U(i,j+1))^2); elseif ((i==64) && (j==1))S = S + (U(i,j) - U(i-1,j))/sqrt(epsilon^2 + (U(i,j) - U(i-1,j))^2);S = S + (U(i,j) - U(i,j+1))/sqrt(epsilon^2 + (U(i,j) - U(i,j+1))^2); elseif ((i==1) &&(j==64))S = S + (U(i,j) - U(i+1,j))/sqrt(epsilon^2 + (U(i,j) - U(i+1,j))^2);S = S + (U(i,j) - U(i,j-1))/sqrt(epsilon^2 + (U(i,j) - U(i,j-1))^2); elseif ((i==64) && (j==64))S = S + (U(i,j) - U(i-1,j))/sqrt(epsilon^2 + (U(i,j) - U(i-1,j))^2);S = S + (U(i,j) - U(i,j-1))/sqrt(epsilon^2 + (U(i,j) - U(i,j-1))^2); elseif ((i==1) && (j>1 && j<64))S = S + (U(i,j) - U(i+1,j))/sqrt(epsilon^2 + (U(i,j) - U(i+1,j))^2);S = S + (U(i,j) - U(i,j-1))/sqrt(epsilon^2 + (U(i,j) - U(i,j-1))^2);S = S + (U(i,j) - U(i,j+1))/sqrt(epsilon^2 + (U(i,j) - U(i,j+1))^2); elseif ((i==64) && (j>1 && j<64))S = S + (U(i,j) - U(i-1,j))/sqrt(epsilon^2 + (U(i,j) - U(i-1,j))^2);S = S + (U(i,j) - U(i,j-1))/sqrt(epsilon^2 + (U(i,j) - U(i,j-1))^2);S = S + (U(i,j) - U(i,j+1))/sqrt(epsilon^2 + (U(i,j) - U(i,j+1))^2); elseif ((j==1) && (i>1 && i<64))S = S + (U(i,j) - U(i-1,j))/sqrt(epsilon^2 + (U(i,j) - U(i-1,j))^2);S = S + (U(i,j) - U(i+1,j))/sqrt(epsilon^2 + (U(i,j) - U(i+1,j))^2);S = S + (U(i,j) - U(i,j+1))/sqrt(epsilon^2 + (U(i,j) - U(i,j+1))^2); elseif ((j==64) && (i>1 && i<64))S = S + (U(i,j) - U(i-1,j))/sqrt(epsilon^2 + (U(i,j) - U(i-1,j))^2);S = S + (U(i,j) - U(i+1,j))/sqrt(epsilon^2 + (U(i,j) - U(i+1,j))^2);S = S + (U(i,j) - U(i,j-1))/sqrt(epsilon^2 + (U(i,j) - U(i,j-1))^2); elseif ((j>1 && j<64) && (i>1 && i<64))S = S + (U(i,j) - U(i-1,j))/sqrt(epsilon^2 + (U(i,j) - U(i-1,j))^2);S = S + (U(i,j) - U(i+1,j))/sqrt(epsilon^2 + (U(i,j) - U(i+1,j))^2);S = S + (U(i,j) - U(i,j-1))/sqrt(epsilon^2 + (U(i,j) - U(i,j-1))^2);S = S + (U(i,j) - U(i,j+1))/sqrt(epsilon^2 + (U(i,j) - U(i,j+1))^2); end f(1) = (W(i,j) - U(i,j) - lambda*S - X(i,j))/mass; f(2) = X(i,j); end

然后寫一下RK4的代碼，$\lambda$是平滑參數，我就嘗試了一下8、12、16，輸入$x_0$是初始值（64*64的矩陣），$h$是步長。

function [y, DY] = RK4(x0,h) x = 0; y = x0; DY = zeros(64,64); delta = 1; k_1 = zeros(64,64,2); k_2 = zeros(64,64,2); k_3 = zeros(64,64,2); k_4 = zeros(64,64,2); N = 0; while (norm(delta)>1e-2 && N<360) % 360 is maximal number of iteration for i = 1:64for j = 1:64k_1(i,j,:) = RHS(x,y,DY,i,j);endendfor i = 1:64for j = 1:64 k_2(i,j,:) = RHS(x+0.5*h,y+0.5*h*squeeze(k_1(:,:,1)),DY+0.5*h*squeeze(k_1(:,:,2)),i,j);endendfor i = 1:64for j = 1:64k_3(i,j,:) = RHS(x+0.5*h,y+0.5*h*squeeze(k_2(:,:,1)),DY+0.5*h*squeeze(k_2(:,:,2)),i,j);endendfor i = 1:64for j = 1:64k_4(i,j,:) = RHS(x+h,y+h*squeeze(k_3(:,:,1)),DY+h*squeeze(k_3(:,:,2)),i,j);endenddelta = (1/6)*(squeeze(k_1(:,:,1))+2*squeeze(k_2(:,:,1))+2*squeeze(k_3(:,:,1))+squeeze(k_4(:,:,1)))*h;y = y + delta; DY = DY + (1/6)*(squeeze(k_1(:,:,2))+2*squeeze(k_2(:,:,2))+2*squeeze(k_3(:,:,2))+squeeze(k_4(:,:,2)))*h; x = x+h;N = N+1; end end

如果$\lambda$取8，輸出是長這樣的

個人感覺幾乎沒啥區別。這張是$\lambda$取12的輸出

這個要稍微白一點了。再看看$\lambda =16$的

這個好像更白一點。這個模型做出來效果其實很一般，要想結果好一點可以再調一下步長，迭代次數上限，平滑參數，停止準則之類，但我就應付個作業，就寫個大概的思想。

總結

以上是生活随笔為你收集整理的UA MATH575B 数值分析下III 图像恢复的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。