CVPR2017的一篇論文
Learning Deep CNN Denoiser Prior for Image Restoration：
一般的，image restoration(IR)任務旨在從觀察的退化變量$y$(退化模型，如式子1)中，恢復潛在的干凈圖像$x$
$y \text{} =\text{}\textbf{H}x\text{}+\text{}v $
where $\textbf{H}$denotes 退化矩陣，$\textbf{v}$denotes 加性高斯白噪聲(additive white Gaussian noise) with 標准差$\sigma$
指定不同的退化矩陣$\textbf{H}$，對應着不同的IR任務：
-- 當$\textbf{H}$是一個恆等矩陣，IR任務對應着圖像去噪(image denoising)
-- 當$\textbf{H}$是一個模糊算子(blurring operator)，IR任務對應着圖像去模糊(image deblurring)
-- 當$\textbf{H}$是一個模糊和下采樣的復合算子(composite operator of blurring and down-sampling)，IR任務對應着圖像超分辨率(image super-resolution)
IR 是一個病態逆問題(ill-posed inverse problem)，先驗(prior)被叫做正則化項(regularization)，需要被采取去約束解決空間. 從貝葉斯的觀點， 需要被推導出的潛在的干凈圖像$\hat{x}$能夠通過解決一個MAP(maximum A posteriori )問題得到:
$\hat{x}\text{}=\text{}$ argmax$_x \text{}$log$p(\textbf{y}|\textbf{x})\text{}+\text{}$log$p(\textbf{x})\text{}\text{}\text{}(1)$
其中log$p(\textbf{y}|\textbf{x})$代表着$y$的log-似然，log$p(x)$是$x$的先驗，更加正式的，公式一 被重新寫：
$\hat{x} = $argmin$_{x}\frac{1}{2}||\textbf{y} - \textbf{H}x||^{2} + \lambda\Phi(x)\text{}\text{}\text{}(2)$
其中解法是最小化一個能量函數，由三個部分組成 數據保真項=$\frac{1}{2}||y\text{}\text{}-\textbf{H}x||^{2}$，一個正則化項$\Phi(x)$ 和一個平衡(trade-off)參數$\lambda$. 數據保真項確保了這個解法是通過了這個退化過程，正則化項強制執行輸出所需的屬性；
To solve Eqn.(2) ， have two main categories methods:
--model-based optimization methods: aim to directly solve Eqn2 使用一些耗時的迭代傳播(time-comsuming iterative inference)
--discriminative learning methods: 通過優化包含退化干凈(degraded-clean pairs )圖像對的訓練集上的損失函數來學習先驗參數$\Theta$和緊湊(compact)推理,目標被給出：
min$_{\Theta}\text{}\textit{l}(\hat{x},x)$ $s.t.$ $\hat{x}\text{}=\text{}$arg min$_{x}\frac{1}{2}||y - \textbf{H}x||\text{}+\text{}\lambda\Phi(x;\Theta))$ (3)
因為這個推理是被MAP估計所指導，所以我們將這種方法稱作為MAP推理指導的判別式學習方法(MAP inference guided discriminative learning method)；
使用一個預先定義的非線性函數:
$\hat{x}\text{}\text{}f(y,\textbf{H};\Theta)$==》注意這種函數內分號的形式，$\Theta$只是參數，前面那個才是真正的自變量,f(x;θ)是關於x的函數，其中θ是參數，x才是自變量，所以f(x;θ)是關於x的一元函數.
來替換MAP推理，我們可以把朴素的判別學習方法看作是Eqn3的一般情況。
--model-based optimization methods: 通過指定退化矩陣$\textbf{H}$來靈活的處理不同的IR任務
--discriminative learning methods: 通過使用確定的退化矩陣訓練數據來學習模型
SO====》由以上兩點可以得到，判別式學習方法常常限定在一些特定的任務如 MLP SRCNN DCNN；
圖像降噪認為model-based method 能夠處理不同噪聲水平的圖像,如BM3D WNNM， 然而判別學習方法只能夠針對不同的噪聲水平訓練不同的模型。
使用判別式學習方法的優缺點：
缺點： 靈活性較差
優點： 快速的測試速度； 由於聯合優化(joint optimization)和端到端的訓練(end-to-end training)，具有前景性的表現；
使用model-based optimization 方法的優缺點：
缺點：時間消耗time-consuming with 復雜的先驗 sophisticated prior 為了好的表現.
SO ===》 非常有吸引力的是利用他們各自的優點去研究他們的集成
Method
在變量分割技術的輔助下，如交替方向乘子法(alternating direction method of multiplier，ADMM)和半二次分割法(half quadratic splitting)，使得有可能獨立的處理 保真項(fidelity term)和正則化項(regularization term).
而且，正則化項僅僅對應於一個去噪子問題！！！所以，這使得能夠將任何的判別式降噪器集成到model-based 優化方法中。
--不采用學習一個MAP推導指導的判別模型，而是采取一個朴素的CNN去學習降噪器；
-- 已學習的CNN降噪器作為一個模塊嵌入到model-based 優化方法中，去解決inverse 問題；
Image Restoration with Denoiser prior
--之前有借助ADMM變量分離技術的即插即用的先驗架構；
--之前有提出使用HQS技術進行圖像的降噪和去模糊的；
基於以上一些降噪先驗的related work, 我們可以看出，降噪器先驗可以通過多種方法插入迭代方案(iterative scheme);
迭代方案涉一個保真項相關子問題和一個降噪子問題；
Method：
--可以幫助噪聲先驗(denoiser prior),which作為基於模型的最優化方法的其中一個模塊來解決這些逆問題(e.g., deblurring).
--噪聲先驗通過判別式學習方法(discriminative learning method)獲得;
So, 結合上面兩點，==》通過CNN訓練一個噪聲器，加入到基於模型的最優化方法來解決其他的逆問題;
在變量分離技術的幫助下，我們可以同時使用兩種方法的各自優點;
變量分割技術(variable spitting techniques):
變量分離技術(variable splitting technique)，如ADMM(alternating direction method of multipliers )，HQS(half quadratic splitting)方法，使得可以分別處理保真項(fidelity term)和正則項(regularization term)，其中正則項僅對應於去噪的子問題。因此，可以在基於模型的優化方法中使用discriminative denoisers,本文的目標在於訓練一系列快速高效的discriminative denoisers，並把它們用於基於模型優化的方法中，解決求逆問題。不使用MAP相關方法，而是使用CNN學習denoisers。
(也可以理解為基於模型的方法一般需要反復迭代去解這個公式，而基於判別學習的方法則通過損失函數去學習先驗參數。這里可以將兩者進行結合，正則項可以對應於一個去噪的子問題，這個子問題可以通過判別式學習的去噪器去獲得，從而帶來圖像先驗，使得基於模型的方法可以快速工作)
貢獻：
-訓練出一系列CNN denoisers。使用變量分離技術，強大的denoisers可以為基於模型的優化方法帶來圖像先驗。
-學習到的CNN denoisers被作為一個模塊部分插入基於模型的優化方法中，解決其他的求逆問題。
半二次方分裂 Half Quadratic Splitting (HQS)
$\hat{x}=\left. arg \text{ }min \right|_{x} \frac{1}{2}|| y-H x||^{2}+\lambda \Phi(x) \text{ }\text{ }\text{ }\text{ }\text{ }\text{ }(2)$
引入輔助變量$z, z = x$,HQS嘗試最小化下面的成本函數:
$L_{\mu}(x,z)=\frac{1}{2}|| y-H x||^{2}+\lambda\Phi(z)+\frac{\mu}{2}||z-x||^{2}\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }(5)$
$\mu$ 懲罰參數，在接下來的迭代求解過程中以非下降階(non-descending order)的形式變化;
等式(5)可以被下面兩個迭代的式子所解決：變量分割技術，
$x_{k+1}=\left. arg\text{ }min \right|_{x}|| y-H x||^{2}+\mu||x-z_k||^{2}\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }(6a)$
$z_{k+1}=\left. arg\text{ } min \right|_{z}\frac{\mu}{2} ||z-x_{k+1}||+\lambda\Phi(z)\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }(6b)$
可以看到保真項與正則化項被分開到兩個子問題中
等式(6a)保真項在二次正則化最小二乘問題，有很多針對不同的退化矩陣的快速解法,最簡單的解法是
$x_{k+1}=(H^{T}H+\mu I)^{-1}(H^{T}y+\mu z_{k}) \text{ }\text{ }\text{ }\text{ }\text{ }\text{ }(7)$
正則化項涉及在6(a)中，可以重寫為(8)
$z_{k+1}=\left. arg\text{ } min \right|_{z} \frac{1}{2(\sqrt{\lambda / \mu})^2} ||x_{k+1}-z||^2+\Phi(z)\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }(8)$
通過貝葉斯概率公式，等式(8)可以看做是對應於一個去噪任務，噪聲水平為$\sqrt{\lambda / \mu}$，所以可以通過去噪器實現求出$z_{k+1}$.
以噪聲水平$\sqrt{\lambda / \mu}$高斯去噪器的去噪圖像$x_{k+1}$.去噪器可以作為(2)的模塊，為了強調這個，重寫(8)
$z_{k+1}=Denoiser(x_{k+1},\sqrt{\lambda / \mu})\text{ }\text{ }\text{ }\text{ }\text{ }\text{ }(9)$
值得注意的是圖像先驗$\Phi$可以間接被去噪先驗替代，這種解法有一些優點：
-- 他允許使用各種灰度和彩色降噪器去解決各種inverse 問題；
-- 求解Eqn2時，顯式圖像先驗$\Phi(\cdot)$是未知的;
-- 利用多個互補(complementary)的去噪器，利用不同的圖像先驗，可以共同解決一個特定的問題;
圖像恢復的MAP推理公式：
$\hat{x}\text{}=\text{}$arg min$_{x}\frac{1}{2}||\textbf{y}\text{}-\text{}\textbf{H}x||^{2}\text{}+\text{}\lambda\Phi(x)$
正則化項$\Phi(x)$對應恢復的表現扮演了至關重要的角色：
$\textbf{z}_{k+1}\text{}=\text{}Denoiser(\textbf{x}_{k+1},\sqrt{\lambda/\mu})$
然后介紹現在的降噪先驗只要采取model-based 優化方法去解決inverse problem，包括：
-- total variation(TV)法 ==》 常常制造watercolor-like 鬼影、偽影
-- 高斯混合模型(Gaussian mixture model，GMM)
-- K-SVD ==》高計算消耗
-- 非局部均值(Non-local means) ==》如果圖像不具有自相似屬性，會過度平滑不規則的結構
-- BM3D==》 如果圖像不具有自相似屬性，會過度平滑不規則的結構
圖像的顏色先驗是一個十分重要的考慮因素，因為圖像大多數圖像是RGB格式。
而由於不同圖像通道之間的相關性，聯合處理圖像的不同通道常常會產生更好的表現比獨立處理每個顏色通道。
許多工作都只對灰度圖像進行建模，而對於彩色圖像的建模較少；
作者指出CBM3D 因為聯合處理了RGB通道，收獲了不錯的效果，同時作者提出可以使用判別學習方法去自動化的揭示潛在的彩色圖像先驗，而不是依靠手工設計的pipeline；
CNN降噪先驗具有速度、表現、判別彩色圖像建模的優勢，同時CNN去學習判別式降噪器(discriminative denoiser)有一些原因：
-- CNN 的前向傳播由於GPU的存在而並行計算
-- CNN 表現出了強大的先驗建模能力with deep architecture
-- CNN利用外部先驗，作為了BM3D為代表的內部先驗的補充
-- 利用判別式學習的優勢
模型結構
CNN denoiser 如上圖所示，網絡包含七層，其中第一層是"擴張卷積(擴張指標為1，感知域還是3*3)+RELU"，2-6層為“擴張卷積(擴張指標分別為2 3 4 3 2)+BN+RELU”， 最后一層為“擴展卷積(1)，相當於正常的卷積運算，且每個中間層的特征圖的數量都為64
擴張卷積filter and 增大的感知域
-- 在圖像降噪中，上下文的信息能夠促進毀壞像素的重建；
-- 為了捕獲上下文的信息，通過前向的卷積操作去增大感知域是一個成功的方法；
有兩種基本的增大感知域的方法：
-- 一個是增大filter size 弊端：會引入更多的參數，增大了計算負擔
-- 一個是增大模型的深度
使用擴張卷積去獲取filter size 和模型深度的平衡， 在保持 3*3 filter 的基礎上 增大感知域， 整個7層網絡實現了 33*33 的感知域， 相當於16層的3*3普通卷積；
其中擴張卷積的filter size 和 擴張指標s 之間的關系為： size = (2s+1)*(2s+1)
使用BN和殘差學習加快訓練
對於高斯降噪問題，結合BN和殘差學習是十分有幫助的，他們都能夠互相的獲益(在他的論文中有講Residual learning of deep CNN for image denoising.） 對於模型的遷移也有用；
殘差學習方式，就是模型的目標不是直接學習產生降噪的圖片，而是學習噪聲即殘差，即輸入的帶噪聲的圖片和干凈圖片的差。
使用小尺寸的圖像作為訓練集去避免邊緣偽影
-- 由於CNN的特點，如果沒有合適的處理，CNN的降噪圖片將會產生邊緣偽影；
-- 對稱pandding 和 zero padding是兩種解決這個問題的方法
-- 對於擴張指數為4的操作，在邊緣pads 4 zeros,那其他的擴張指數呢？
-- 經驗主義的使用了小尺寸的訓練樣本 去避免邊緣偽影，原因包括：
-- 將大尺寸的圖像crop 成小尺寸的patches,有利於CNN去看到更多的邊沿信息，比如將70*70的patches crop成四個非重疊的35*35的patches,邊緣信息被擴大了；
-- 就patch的大小可以作對比試驗進行驗證；
-- 當訓練的patch尺寸小於感知域后，這個性能會下降；
學習實際的降噪器模型with 小間隔的噪聲水平
-- 想要得到精確的子問題的解是非常困難且time-consuming的去優化的，使用不精確但是快速的子問題的解能夠加快收斂(兩篇文獻：The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices.和 From learning models of natural im age patches to whole image restoration.)
-- 所以 沒有必要去學習很多判別式降噪模型for 每個噪聲水平。
-- 盡管$\textbf{z}_{k+1}\text{}=\text{}Denoiser(\textbf{x}_{k+1},\sqrt{\lambda/\mu})$ 是一個降噪器，但他與傳統的高斯降噪有着不同的目標。
-- 傳統的高斯降噪是恢復出潛在的干凈圖像，無論要去噪的圖像的噪聲類型和噪聲水平如何，這里的去噪器都會發揮自己的作用。也就是說，不管這個圖像有沒有噪聲，都會發揮作用！
-- 所以一個理想的判別降噪器應該使用當前的噪聲水平進行訓練：
訓練了一系列的噪聲水平在0-50同時獨立的以2為間隔的模型，產生了25個模型為圖像的先驗進行建模；迭代方案的存在，使他恢復足以滿足。
實驗
圖像降噪：
-- 將每個圖片 crop 成了35*35的patches,因為使用殘差學習的方式,損失函數：
$\textit{l}(\Theta)\text{}=\text{}\frac{1}{2N}\sum_{i=1}^{N}||f(y_{i};\Theta)\text{}-\text{}(y_{i}\text{}-\text{}x_{i})||^{2}$
-- 訓練結束的標志 訓練損失在五個連續的epoch固定
-- 使用了旋轉翻轉等數據擴充技巧；
-- 從不同方法的PSNR 進行了對比， 分別書灰度圖和彩色圖
-- 從不同方法的運行時間進行了對比
圖像去模糊：
模糊核的選擇：
-- 一個常見的模糊的高斯模糊核，標准差為1.6， 來自論文(Understanding and evaluating blind deconvolution algorithms)的前兩個的真實模糊核；
-- 我們只需將顏色去噪器插入到HQS框架中；
-- 在公式6中：
-- 兩個參數中，$\lambda$ 與$\sigma^{2}$相聯系同時在迭代中保持固定，其中$\mu$控制着降噪器的噪聲水平；
公式6(a)的快速解法：
$x_{k+1}=(H^{T}H+\mu I)^{-1}(H^{T}y+\mu z_{k}) \text{ }\text{ }\text{ }\text{ }\text{ }\text{ }(7)$
-- 由於hqs框架是基於去噪的，因此我們將每次迭代中去噪的噪聲級隱式的確定為μ。
-- 在我們的實驗設置中，根據噪聲水平，它以指數形式從49衰減到[1，15]中的值。
-- 實驗代碼核心迭代部分
%模擬仿真的模糊噪聲圖片
k = fspecial('gaussian', 25, 1.6);%相當於是模糊算子
y = imfilter(im2double(x), k, 'circular', 'conv') + sigma*randn(size(x));
%denominator 表示是H^{T}H, H模糊算子與他的共軛矩陣相乘
V = psf2otf(k,[w,h]);%psf2otf(PSF) 其作用是將一個空間點擴散函數轉換為頻譜面的光學傳遞函數，執行的也是對PSF的FFT變換，變為了頻域；
denominator = abs(V).^2;%y = abs(3+4i) y=5
%H^{T}*y H的共軛矩陣乘以輸入y
upperleft = conj(V).*fft2(y); % conj(V) V=5-2i ==> 5+2i
%訓練的特點噪聲水平\sigma的一組denoiser，其\sigma以指數形式的衰減的不同模型；
modelSigmaS = logspace(log10(modelSigma1),log10(modelSigma2),totalIter);
%\lamba與噪聲水平\sigma^2有關且在迭代中保持固定，\mu控制着降噪器的噪聲水平，%rho隱式的表示\mu，見下面的計算方法與降噪器的噪聲水平和噪聲水平有關
rho = lamda*255^2/(modelSigmaS(itern)^2);%[243.742,267.099,292.694,320.742,351.476]
upperleft =
for itern = 1:totalIter
%%% step 1
rho = lamda*255^2/(modelSigmaS(itern)^2);
z = real(ifft2((upperleft + rho*fft2(z))./(denominator + rho)));
if ns(itern+1)~=ns(itern)
[net] = loadmodel(modelSigmaS(itern),CNNdenoiser);
net = vl_simplenn_tidy(net);
if useGPU
net = vl_simplenn_move(net, 'gpu');
end
end
%%% step 2
res = vl_simplenn(net, z,[],[],'conserveMemory',true,'mode','test');
residual = res(end).x;
z = z - residual;
end
真實模糊圖像的測試
設置了兩個重要的估計圖像噪聲水平 和 降噪器去噪水平的 超參數：
% There are two important parameters to tune:
% (1) image noise level of blurred image: Isigma and
% (2) noise level of the last denoiser: Msigma.
%使用了 別人的方法 產生的模糊核 圖像作為先驗的條件
%% read blurred image and its estimated kernel
% blurred image
Iname = 'im01_ker01';
y = im2single(imread(fullfile(folderTestCur,[Iname,'.png'])));
% estimated kernel
%k = imread(fullfile(folderTestCur,[Iname,'_kernel.png']));
k = imread(fullfile(folderTestCur,[Iname,'_out_kernel.png']));
if size(k,3)==3
k = rgb2gray(k);
end
k = im2single(k);
k = k./(sum(k(:))); %歸一化
%比較重要的部分是邊緣的處理
%% handle boundary
boundary_handle = 'case2';
switch boundary_handle
case {'case1'} % option (1), edgetaper to better handle circular boundary conditions, (matlab2015b)
% k(k==0) = 1e-10; % uncomment this for matlab 2016--2018?
ks = floor((size(k) - 1)/2);
y = padarray(y, ks, 'replicate', 'both');
for a=1:4
y = edgetaper(y, k);
end
case {'case2'} % option (2)
H = size(y,1); W = size(y,2);
y = wrap_boundary_liu(y, opt_fft_size([H W]+size(k)-1));
end
需要明白為什么需要進行邊緣條件的變化，是讓模糊核完全重合和圖像，所以需要用模糊核 進行 一些 padding 操作
降噪的代碼，論文並沒有提供在實際環境中應用的代碼；
% clear; clc;
addpath('utilities');
imageSets = {'CBSD68','Val20'}; %%% testing dataset
folderTest = 'testsets';
folderModel = 'models';
folderResult = 'results';
taskTestCur = 'Denoising';
if ~exist(folderResult,'file')
mkdir(folderResult);
end
setTestCur = imageSets{1};
imageSigmaS = [5,10,15,25,35,50];
modelSigmaS = [5,10,15,25,35,50];
showResult = 1;
saveResult = 0;
useGPU = 1;
pauseTime = 1;
%%% folder to store results
folderResultCur = fullfile(folderResult, [taskTestCur,'_',setTestCur]);
if ~exist(folderResultCur,'file')
mkdir(folderResultCur);
end
%%% read images
ext = {'*.jpg','*.png','*.bmp'};
filePaths = [];
folderTestCur = fullfile(folderTest,setTestCur);
for i = 1 : length(ext)
filePaths = cat(1,filePaths, dir(fullfile(folderTestCur,ext{i})));
end
%%% PSNR and SSIM
PSNRs = zeros(length(modelSigmaS),length(filePaths));
SSIMs = zeros(length(modelSigmaS),length(filePaths));
load(fullfile(folderModel,'modelcolor.mat'));
for i = 1:length(modelSigmaS)
disp([i,length(modelSigmaS)]);
net = loadmodel(modelSigmaS(i),CNNdenoiser);
net = vl_simplenn_tidy(net);
% for i = 1:size(net.layers,2)
% net.layers{i}.precious = 1;
% end
%%% move to gpu
if useGPU
net = vl_simplenn_move(net, 'gpu');
end
for j = 1:length(filePaths)
%%% read images
label = imread(fullfile(folderTestCur,filePaths(j).name));
[~,imageName,extCur] = fileparts(filePaths(j).name);
label = im2double(label);
randn('seed',0);
input = single(label + imageSigmaS(i)/255*randn(size(label)));
%%% convert to GPU
if useGPU
input = gpuArray(input);
end
res = vl_simplenn(net,input,[],[],'conserveMemory',true,'mode','test');
output = input - res(end).x;
%%% convert to CPU
if useGPU
output = gather(output);
input = gather(input);
end
%%% calculate PSNR and SSIM
[PSNRCur, SSIMCur] = Cal_PSNRSSIM(im2uint8(label),im2uint8(output),0,0);
if showResult
imshow(cat(2,im2uint8(label),im2uint8(input),im2uint8(output)));
title([filePaths(j).name,' ',num2str(PSNRCur,'%2.2f'),'dB',' ',num2str(SSIMCur,'%2.4f')])
drawnow;
if saveResult
imwrite(im2uint8(output),fullfile(folderResultCur,[imageName,'_',num2str(imageSigmaS(i)),'_',num2str(modelSigmaS(i)),'_',num2str(PSNRCur,'%2.2f'),'.png']));
end
pause(pauseTime)
end
PSNRs(i,j) = PSNRCur;
SSIMs(i,j) = SSIMCur;
end
end
%%% save PSNR and SSIM metrics
save(fullfile(folderResultCur,['PSNR_',taskTestCur,'_',setTestCur,'.mat']),'PSNRs')
save(fullfile(folderResultCur,['SSIM_',taskTestCur,'_',setTestCur,'.mat']),'SSIMs')
disp([mean(PSNRs,2),mean(SSIMs,2)]);