7、图像复原与重建

作者: sumpig | 来源:发表于2019-03-28 10:47 被阅读0次

图像复原试图利用退化现象的某种先验知识来复原一幅退化的图像，大部分是客观处理。例如通过去模糊函数去除图像模糊则被认为是一种图像复原技术。

1、图像复原模型

退化过程被建模为一个退化函数和一个加性噪声项，退化后的图像为：

$g(x,y)=H[f(x,y)]+\eta(x,y)$

空间域退化公式：
$g(x,y)=h(x,y)★f(x,y)+\eta(x,y)$

频率域退化公式：
$G(u,v)=H(u,v)F(u,v)+N(u,v)$

F(u,v) 光传递函数（OTF）
h(x,y) 点扩散函数（PSF）

2、噪声模型

使用函数 imnoise 对图像添加噪声

通过函数 imnoise 用噪声污染一幅图像，语法为

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 的均匀分布随机噪声

g=imnoise(f, 'poisson')
泊松噪声

使用规定分布生成空间随机噪声

通常，我们需要生成函数 imnoise 所不能生成的噪声类型和参数。空间噪声值是随机数，可采用一些简单的概率论规则来生成分布类型中我们感兴趣的随机数。

下面这个自定义函数可以生成噪声模式本身：

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); % 对数正太
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;

% --------------------- 指数 ------------------- %
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

示例

下面这个语句生成 100000 个元素的列向量 r ，每个元素都是一个随机数，这些随机数服从高斯分布，均值为0，标准差为 1。

r = imnoise2('gaussian', 100000, 1, 0, 1);

可以使用 hist 得到 r 的直方图

hist(r, 50);

untitled.jpg

周期噪声

周期噪声通常源于电气电机干扰，通常是通过在频率域滤波来处理。
我们的周期噪声模型是一个离散的二维正弦波：

$r(x, y) =Asin[2\pi u_0(x+B_x)/M + 2\pi v_0(y+B_y)/N]$

下面的函数接受任意数量的脉冲位置（频率坐标），每个脉冲位置都有自己的增幅、频率和相移参数，该函数还会输出各个正弦波之和的傅里叶变换及相应谱。

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, [])

估计噪声参数（略）

function [p, npix] = histroi(f, c, r)
B = roipoly(f, c, r);
p = imhist(f(B));
if nargout > 1
    npix = sum(B(:));
end

3、仅有噪声的复原 -- 空间滤波

如果退化仅是噪声，那么模型变为

$g(x, y)=f(x,y)+\eta(x,y)$

这种情况下，所选的降噪方法是空间滤波。

空间噪声滤波器

下面的自定义函数可用相应的滤波器执行空间滤波。

function f = spfilt(g, type, varargin)
switch type
case 'amean' % 算术均值
    w = fspecial('average', [m n]);
    f = imfilter(g, w, 'replicate');
case 'gmean' % 几个均值
    f = gmean(g, m, n);
case 'hmean' % 调和均值
    f = harmean(g, m, n);
case 'chmean' % 逆调和均值
    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' % 中点值
    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

自适应空间滤波器（略）

f = adpmedian(g, Smax)

4、退化函数建模（略）

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、直接逆滤波（略）

$\hat F(u,v)=\frac{G(u, v)}{H(u, v)}$

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