UA MATH575B 數(shù)值分析下III 圖像恢復(fù)

• 圖像去噪的優(yōu)化模型
• 算法實(shí)現(xiàn)

圖像去噪的優(yōu)化模型

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

算法實(shí)現(xiàn)

為什么要用這個(gè)枯燥的例子呢？因?yàn)檫@是我的作業(yè)【吐血.jpg】下面這張灰度圖是我要用來(lái)去噪的，不難看出它居然是個(gè)A！

這張圖可以用一個(gè)矩陣來(lái)表示，對(duì)應(yīng)上面的記號(hào)，用$W$來(lái)表示，大概長(zhǎng)下面這樣，每個(gè)值代表那個(gè)像素點(diǎn)的灰度，取值是0-255，每個(gè)值占8 bits。上面的原圖就是用這個(gè)矩陣生成的，在matlab里可以用這個(gè)來(lái)生成

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個(gè)像素點(diǎn)，需要重構(gòu)的圖像$U$是一個(gè)有4096個(gè)元素的矩陣。在這個(gè)ODE-based梯度下降的計(jì)算模型中
$$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}}$$
參考老師給的參數(shù)，$m=0.2,?=0.1$。先把這個(gè)兩階的ODE用二元一階ODE表示一下，假設(shè)$X_v=\frac{d{U}_v}{dt}$，則這個(gè)方程可以寫(xiě)成
$$\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}$$
這個(gè)看起來(lái)是二元一階ODE的東西其實(shí)有8192個(gè)未知數(shù)。現(xiàn)在把這個(gè)圖像讀進(jìn)Matlab，再定義一些參數(shù)

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

定義一下ODE右邊那個(gè)函數(shù)，RK4會(huì)用到，要注意四個(gè)角和四條邊上的像素相鄰像素點(diǎn)分別只有2和3個(gè)。

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

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

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，輸出是長(zhǎng)這樣的

個(gè)人感覺(jué)幾乎沒(méi)啥區(qū)別。這張是$\lambda$取12的輸出

這個(gè)要稍微白一點(diǎn)了。再看看$\lambda =16$的

這個(gè)好像更白一點(diǎn)。這個(gè)模型做出來(lái)效果其實(shí)很一般，要想結(jié)果好一點(diǎn)可以再調(diào)一下步長(zhǎng)，迭代次數(shù)上限，平滑參數(shù)，停止準(zhǔn)則之類，但我就應(yīng)付個(gè)作業(yè)，就寫(xiě)個(gè)大概的思想。

總結(jié)

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

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。