图像复原与重建

一、概念

1.图像恢复
质量变坏，模糊、失真、有噪声。
2.图像复原
使退化的图像尽量恢复本来面目，是图像退化的逆过程进行处理。
3.图像复原的过程
重点是找到图像退化的原因

噪声模型

使用函数 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);

周期噪声

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

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