圖像復(fù)原試圖利用退化現(xiàn)象的某種先驗(yàn)知識(shí)來復(fù)原一幅退化的圖像赖临，大部分是客觀處理胞锰。例如通過去模糊函數(shù)去除圖像模糊則被認(rèn)為是一種圖像復(fù)原技術(shù)。
1兢榨、圖像復(fù)原模型
退化過程被建模為一個(gè)退化函數(shù)和一個(gè)加性噪聲項(xiàng)嗅榕，退化后的圖像為：
空間域退化公式：
頻率域退化公式：
F(u,v) 光傳遞函數(shù)（OTF）
h(x,y) 點(diǎn)擴(kuò)散函數(shù)（PSF）
2、噪聲模型
使用函數(shù) imnoise 對(duì)圖像添加噪聲
通過函數(shù) imnoise 用噪聲污染一幅圖像吵聪，語(yǔ)法為
g=imnoise(f, type, parameters)
% type: gaussian, localvar, salt & pepper, speckle, poisson
- g=imnoise(f, 'gaussian', m, var)
均值為 m凌那、方差為 var 的高斯噪聲- g=imnoise(f, 'localvar', V)
均值為 0、局部方差為 V 的高斯噪聲- g=imnoise(f, 'salt & pepper', d)
噪聲密度為 d 的椒鹽噪聲- g=imnoise(f, 'speckle', var)
均值為 0吟逝、方差為 var 的均勻分布隨機(jī)噪聲- g=imnoise(f, 'poisson')
泊松噪聲
使用規(guī)定分布生成空間隨機(jī)噪聲
通常帽蝶，我們需要生成函數(shù) imnoise 所不能生成的噪聲類型和參數(shù)】樵埽空間噪聲值是隨機(jī)數(shù)励稳，可采用一些簡(jiǎn)單的概率論規(guī)則來生成分布類型中我們感興趣的隨機(jī)數(shù)。
下面這個(gè)自定義函數(shù)可以生成噪聲模式本身：
function R = imnoise2(type, varargin)
[M, N, a, b] = setDefaults(type, varargin{:});
switch lower(type)
case 'uniform'
R = a + (b-a)*rand(M, N); % 均勻
case 'gaussian'
R = a + b*randn(M, N); % 高斯
case 'salt & pepper'
R = saltpepper(M, N, a, b);
case 'lognormal'
R = exp(b*randn(M, N) + a); % 對(duì)數(shù)正太
case 'rayleigh'
R = a + (-b*log(1 - rand(M, N))).^0.5; % 瑞利
case 'exponential'
R = exponential(M, N, a);
case 'erlang'
R = erlang(M, N, a, b);
otherwise
error('Unknow distribution type.')
end
% ------------------- 椒鹽 -------------------- %
function R = saltpepper(M, N, a, b)
if (a + b) > 1
error('The sum Pa + Pb must not exceed 1.')
end
R(1:M, 1:N) = 0.5;
X = rand(M, N);
R(X <= a) = 0;
u = a + b;
R(X > a & X <= u) = 1;
% --------------------- 指數(shù) ------------------- %
function R = exponential(M, N, a)
if a <= 0
error('Parameter a must be positive for exponential type.')
end
k = -1 / a;
R = k * log(1-rand(M, N));
% -------------------- 愛爾蘭 -------------------- %
function R = erlang(M, N, a, b)
if (b ~= round(b) || b <= 0)
error('Param b must be a positive integer for Erlang.')
end
k = -1 / a;
R = zeros(M, N);
for j = 1:b
R = R + k*log(1 - rand(M, N));
end
% --------------------------------------- %
function varargout = setDefaults(type, varargin)
varargout = varargin;
P = numel(varargin);
if P < 4
varargout{4} = 1;
end
if P < 3
varargout{3} = 0;
end
if P < 2
varargout{2} = 1;
end
if P < 1
varargout{1} = 1;
end
if (P <= 2)
switch type
case 'salt & pepper'
varargout{3} = 0.05;
varargout{4} = 0.05;
case 'lognormal'
varargout{3} = 1;
varargout{4} = 0.25;
case 'exponential'
varargout{3} = 1;
case 'erlang'
varargout{3} = 2;
varargout{4} = 5;
end
end
示例
下面這個(gè)語(yǔ)句生成 100000 個(gè)元素的列向量 r 局蚀，每個(gè)元素都是一個(gè)隨機(jī)數(shù)麦锯，這些隨機(jī)數(shù)服從高斯分布，均值為0琅绅， 標(biāo)準(zhǔn)差為 1。
r = imnoise2('gaussian', 100000, 1, 0, 1);
可以使用 hist 得到 r 的直方圖
hist(r, 50);
untitled.jpg
周期噪聲
周期噪聲通常源于電氣電機(jī)干擾鹅巍，通常是通過在頻率域?yàn)V波來處理千扶。
我們的周期噪聲模型是一個(gè)離散的二維正弦波：
下面的函數(shù)接受任意數(shù)量的脈沖位置（頻率坐標(biāo)），每個(gè)脈沖位置都有自己的增幅骆捧、頻率和相移參數(shù)澎羞，該函數(shù)還會(huì)輸出各個(gè)正弦波之和的傅里葉變換及相應(yīng)譜。
function [r, R, S] = imnoise3(M, N, C, A, B)
K = size(C, 1);
if nargin < 4
A = ones(1, K);
end
if nargin < 5
B = zeros(K, 2);
end
R = zeros(M, N);
for j = 1:K
u1 = floor(M/2) + 1 - C(j, 1);
v1 = floor(N/2) + 1 - C(j, 2);
R(u1, v1) = i*M*n*(A(j).2) * exp(-i*2*pi*(C(j, 1)*B(j, 1)/M ...
+ C(j, 2)*B(j, 2)/N));
u2 = floor(M/2) + 1 + C(j, 1);
v2 = floor(N/2) + 1 + C(j, 2);
R(u2, v2) = -i*M*n*(A(j).2) * exp(i*2*pi*(C(j, 1)*B(j, 1)/M ...
+ C(j, 2)*B(j, 2)/N));
end
S = abs(R);
r = real(ifft2(ifftshift(R)));
示例
C = [0 64; 0 128; 32 32; 64 0; 128 0; -32 32];
[r, R, S] = imnoise3(512, 512, C);
imshow(S, [])
figure, imshow(r, [])
估計(jì)噪聲參數(shù)（略）
function [p, npix] = histroi(f, c, r)
B = roipoly(f, c, r);
p = imhist(f(B));
if nargout > 1
npix = sum(B(:));
end
3敛苇、僅有噪聲的復(fù)原 -- 空間濾波
如果退化僅是噪聲妆绞，那么模型變?yōu)?/p>
這種情況下顺呕，所選的降噪方法是空間濾波。
空間噪聲濾波器
下面的自定義函數(shù)可用相應(yīng)的濾波器執(zhí)行空間濾波括饶。
function f = spfilt(g, type, varargin)
switch type
case 'amean' % 算術(shù)均值
w = fspecial('average', [m n]);
f = imfilter(g, w, 'replicate');
case 'gmean' % 幾個(gè)均值
f = gmean(g, m, n);
case 'hmean' % 調(diào)和均值
f = harmean(g, m, n);
case 'chmean' % 逆調(diào)和均值
f = charmean(g, m, n, Q);
case 'mdian' % 中值
f = medfilt2(g, [m n], 'symmetric');
case 'max' % 最大值
f = imdilate(g, ones(m, n));
case 'min' % 最小值
f = imerode(g, ones(m, n));
case 'midpoint' % 中點(diǎn)值
f1 = ordfilt2(g, 1, ones(m, n), 'symmetric');
f2 = ordfilt2(g, m*n, ones(m, n), 'symmetric');
f = imlicomb(0.5, f1, 0.5, f2);
case 'atrimmed' % 修正的 a 均值
f = alphatrim(g, m, n, d);
otherwise
error('Unknown filter type.')
end
% ------------------------------------ %
function f = gmean(g, m, n)
[g, revertClass] = tofloat(g);
f = exp(imfilter(log(g), ones(m, n), 'replicate')) .^ (1/m/n);
f = revertClass(f);
% ------------------------------------ %
function f = harmean(g, m, n)
[g, revertClass] = tofloat(g);
f = m * n ./ imfilter(1 ./(g+eps), ones(m, n), 'replicate');
f = revertClass(f);
% ------------------------------------ %
function = charmean(g, m, n, q)
[g, revertClass] = tofloat(g);
f = imfilter(g .^ (q+1), ones(m, n), 'replicate');
f = revertClass(f);
% ------------------------------------ %
function f = alphatrim(g, m, n, d)
if (d <= 0) || (d/2 ~= round(d/2))
error('d must be a positive, even interger.')
end
[g, revertClass] = tofloat(g);
f = imfilter(g, ones(m, n), 'symmetric');
for k = 1:d/2
f = f - ordfilt2(g, k, ones(m, n), 'symmetric');
end
for k = (m*n - (d/2) + 1):m*n
f = f - ordfilt2(g, k, ones(m, n), 'symmetric');
end
f = f / (m*n - d);
f = revertClass(f);
% ----------------------------------- %
function [m, n, Q, d] = processInputs(varargin)
m = 3;
n = 3;
Q = 1.5;
d = 2;
if nargin > 0
m = varargin{1};
end
if nargin > 1
n = varargin{2};
end
if nargin > 2
Q = varargin{3};
d = varargin{3};
end
自適應(yīng)空間濾波器（略）
f = adpmedian(g, Smax)
4株茶、退化函數(shù)建模（略）
f = checkerboard(8);
PSF = fspecial('motion', 7, 45);
gb = imfilter(f, PSF, 'circular');
noise = imnoise2('Gaussian', size(f, 1), size(f, 2), 0, sqrt(0.001));
g = gb + noise;
5、直接逆濾波（略）
6图焰、維納濾波（略）
- deconvmnr
7启盛、由投影重建圖像（略）
- radon
g1 = zeros(600, 600);
g1(100:500, 250:350) = 1;
g2 = phantom('Modified Shepp-Logan', 600);
theta = 0:0.5:179.5;
[R1, xp1] = radon(g1, theta);
[R2, xp2] = radon(g2, theta);
R1 = flipud(R1');
R2 = flipud(R2');
imshow(R1, [], 'XData', xp1([1 end]), 'YData', [179.5 0])
- iradon
- fanbeam
