7、图像复原与重建

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

7、图像复原与重建
数字图像处理(五) 图像复原
数字图像处理复习（三）
1. 图像处理、计算机视觉与OpenCV
去雾研究
exp3-图像复原
基于改进GFPGAN的模糊人脸修复系统（源码＆教程）
信号处理(一)
自然图像先验
复原c2018-11-13

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

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

7、图像复原与重建
图像复原试图利用退化现象的某种先验知识来复原一幅退化的图像，大部分是客观处理。例如通过去模糊函数去除图像模糊则被认...
数字图像处理(五) 图像复原
本节主要目的是介绍图像复原一些基本概念，如图像退化/复原过程的模型，图像复原的滤波方法，包括非约束复原(逆滤波...
数字图像处理复习（三）
图像复原基本概念图像复原和图像增强的区别图像增强：是一个主观的过程，从视觉角度，改善图像质量。图像增强被认为...
1. 图像处理、计算机视觉与OpenCV
图像处理（数字图像处理）：图像处理是用计算机对图像进行分析，以达到所需结果的技术，主要包括图像压缩，增强与复原，...
去雾研究
基于物理的模型的方法（图像复原）利用大气散射模型：通过求解图像降质的逆过程来恢复清晰的图像，属于图像复原的范畴2...
exp3-图像复原
图像复原中的均值滤波、统计排序滤波，去除高斯、椒盐噪声去除周期噪声由退化函数进行图像复原
基于改进GFPGAN的模糊人脸修复系统（源码＆教程）
1.研究背景图像是人类获得信息的重要来源，因此图像复原是一项很重要的技术。图像复原指的是修复图像缺失的部分。随着...
信号处理(一)
本篇介绍图像处理离不开采样与重建，本篇就介绍下采样与重建背后的数学逻辑。一维采样采样就是将模拟信号用数字信号...
自然图像先验
自然图像先验在自然图像处理领域里,有很多问题(比如图像去噪、图像去模糊、图像修复、图像重建等)都是反问题 ,即...
复原c2018-11-13
由投影重建图像：滤波反投影、FDK、TFDK三维重建算法理论基础