图像恢复

图像降指与恢复过程

Paste_Image.png

噪声模型

imnoise函数生成噪声

用之前需要先将f规范化
g = imnoise(f, type, parameters)

type:

'gussian' 'salt & pepper' 'motion'

parameters

生成特定分布的空域随机噪声

使用累积分布的逆函数，可以生成任意你想要的原始分布；这里假设w是均匀分布，可以得到z的分布；

Paste_Image.png
Paste_Image.png

CDF是从PDF积分得到的

Paste_Image.png

% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% 模拟噪声
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
% 
close all;
clear all;
r={};
r{1}=imnoise2('gaussian',100000,2,0,1);
r{2}=imnoise2('uniform',100000,2,0,1);
r{3}=imnoise2('rayleigh',100000,2,0,1);

% for i=1:3
%     subplot(2,2,i);imshow(r{i},50);
% end

% figure;
for i=1:numel(r)
   subplot(2,2,i);mesh(r{i}); 
end

figure;
for i=1:numel(r)
   subplot(2,2,i);hist(r{i});
   size(r{i})
end

只存在空域噪声的恢复

与图像增强的技术一致

均值滤波器

Paste_Image.png

统计排序与调和均值滤波器

Paste_Image.png

空域滤波基础函数Spfilt

function f = spfilt(g, type, m, n, parameter)
%SPFILT Performs linear and nonlinear spatial filtering.
%   F = SPFILT(G, TYPE, M, N, PARAMETER) performs spatial filtering
%   of image G using a TYPE filter of size M-by-N. Valid calls to
%   SPFILT are as follows: 
%
%     F = SPFILT(G, 'amean', M, N)       Arithmetic mean filtering.
%     F = SPFILT(G, 'gmean', M, N)       Geometric mean filtering.
%     F = SPFILT(G, 'hmean', M, N)       Harmonic mean filtering.
%     F = SPFILT(G, 'chmean', M, N, Q)   Contraharmonic mean
%                                        filtering of order Q. The
%                                        default is Q = 1.5.
%     F = SPFILT(G, 'median', M, N)      Median filtering.
%     F = SPFILT(G, 'max', M, N)         Max filtering.
%     F = SPFILT(G, 'min', M, N)         Min filtering.
%     F = SPFILT(G, 'midpoint', M, N)    Midpoint filtering.
%     F = SPFILT(G, 'atrimmed', M, N, D) Alpha-trimmed mean filtering.
%                                        Parameter D must be a
%                                        nonnegative even integer;
%                                        its default value is D = 2.
% 
%   The default values when only G and TYPE are input are M = N = 3,
%   Q = 1.5, and D = 2. 

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.6 $  $Date: 2003/10/27 20:07:00 $

% Process inputs.
if nargin == 2
   m = 3; n = 3; Q = 1.5; d = 2;
elseif nargin == 5
   Q = parameter; d = parameter;
elseif nargin == 4
   Q = 1.5; d = 2;
else 
   error('Wrong number of inputs.');
end

% Do the filtering.
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 'median'
   f = medfilt2(g, [m n], 'symmetric');
case 'max'
   f = ordfilt2(g, m*n, ones(m, n), 'symmetric');
case 'min'
   f = ordfilt2(g, 1, ones(m, n), 'symmetric');
case 'midpoint'
   f1 = ordfilt2(g, 1, ones(m, n), 'symmetric');
   f2 = ordfilt2(g, m*n, ones(m, n), 'symmetric');
   f = imlincomb(0.5, f1, 0.5, f2);
case 'atrimmed'
   if (d <= 0) | (d/2 ~= round(d/2))
      error('d must be a positive, even integer.')
   end
   f = alphatrim(g, m, n, d);
otherwise
   error('Unknown filter type.')
end

%-------------------------------------------------------------------%
function f = gmean(g, m, n)
%  Implements a geometric mean filter.
inclass = class(g);
g = im2double(g);
% Disable log(0) warning.
warning off;
f = exp(imfilter(log(g), ones(m, n), 'replicate')).^(1 / m / n);
warning on;
f = changeclass(inclass, f);

%-------------------------------------------------------------------%
function f = harmean(g, m, n)
%  Implements a harmonic mean filter.
inclass = class(g);
g = im2double(g);
f = m * n ./ imfilter(1./(g + eps),ones(m, n), 'replicate');
f = changeclass(inclass, f);

%-------------------------------------------------------------------%
function f = charmean(g, m, n, q)
%  Implements a contraharmonic mean filter.
inclass = class(g);
g = im2double(g);
f = imfilter(g.^(q+1), ones(m, n), 'replicate');
f = f ./ (imfilter(g.^q, ones(m, n), 'replicate') + eps);
f = changeclass(inclass, f);

%-------------------------------------------------------------------%
function f = alphatrim(g, m, n, d)
%  Implements an alpha-trimmed mean filter.
inclass = class(g);
g = im2double(g);
f = imfilter(g, ones(m, n), 'symmetric');
for k = 1:d/2
   f = imsubtract(f, ordfilt2(g, k, ones(m, n), 'symmetric'));
end
for k = (m*n - (d/2) + 1):m*n
   f = imsubtract(f, ordfilt2(g, k, ones(m, n), 'symmetric'));
end
f = f / (m*n - d);
f = changeclass(inclass, f);

逆谐波函数用于消除椒盐噪声

% 逆调和均值滤波的例子

gp  = imread('E:\资料\onedrive\code\test\image\Fig0505(a)(ckt_pepper_only).tif');
gme = spfilt(gp,'gmean',3,3);
gch = spfilt(gp,'chmean',3,3,1.5); % 椒噪声用正的Q;
gar = spfilt(gp,'amean',3,3); 
gme = spfilt(gp,'gmean',3,3);
myImshow(gp,'椒噪');myImshow(gch,'逆调和除噪');myImshow(gar,'均值滤波');%gme
myImshow(gme,'几何均值滤波');
myImshow(0);
figure;
gs  = imread('E:\资料\onedrive\code\test\image\Fig0505(b)(ckt_salt_only).tif');
gch = spfilt(gs,'chmean',3,3,-1.5); % 椒噪声用正的Q;
gar = spfilt(gs,'amean',3,3); 
gme = spfilt(gs,'gmean',3,3);
myImshow(gs,'盐噪');myImshow(gch,'逆调和除噪');myImshow(gar,'均值滤波');
myImshow(gme,'几何均值滤波');
myImshow(0);

untitled.png untitled.png

自适应噪声滤波器

% 中值滤波与自适应中值滤波
gsp = imread('E:\资料\onedrive\code\test\image\Fig0506(a)(ckt_salt_pep_pt25).tif');%Fig0506(a)(ckt_salt_pep_pt25).tif
gm2 = medfilt2(gsp);% 中值滤波
gch1 = spfilt(gsp,'chmean',3,3,1.5); % 椒噪声用正的Q;
gch2 = spfilt(gsp,'chmean',3,3,-1.5); % 盐噪声用负的Q;
gch3 = spfilt(gch1,'chmean',3,3,-1.5); % 除椒后除盐
gadp = adpmedian(gsp,7);
myImshow(gsp,'椒盐噪声');myImshow(gm2,'中值滤波');myImshow(gadp,'自适应中值');
myImshow(gch1,'除椒');myImshow(gch2,'除盐');myImshow(gch3,'除椒后除盐');
myImshow(0);

untitled.png

运动模型的退化函数

运动模型

在快门曝光的瞬间内，取景区与镜头发生了位移，使得原始图片f(x,y)变成了g(x,y)

Paste_Image.png
Paste_Image.png
Paste_Image.png

基础函数

function B = pixeldup(A, m, n)
%PIXELDUP Duplicates pixels of an image in both directions.
%   B = PIXELDUP(A, M, N) duplicates each pixel of A M times in the
%   vertical direction and N times in the horizontal direction.
%   Parameters M and N must be integers.  If N is not included, it
%   defaults to M.

%   Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
%   Digital Image Processing Using MATLAB, Prentice-Hall, 2004
%   $Revision: 1.4 $  $Date: 2003/06/14 16:29:54 $

% Check inputs.
if nargin < 2 
   error('At least two inputs are required.'); 
end
if nargin == 2 
   n = m; 
end

% Generate a vector with elements 1:size(A, 1).
u = 1:size(A, 1);

% Duplicate each element of the vector m times.
m = round(m); % Protect against nonintergers.
u = u(ones(1, m), :);
u = u(:);

% Now repeat for the other direction.
v = 1:size(A, 2);
n = round(n);
v = v(ones(1, n), :);
v = v(:);
B = A(u, v);

u = u(ones(1, m), :);这代码可以推广到一般的形式A(a,b);其中A为矩阵，a、b为向量，则a表示从A中取元素时的所取的行下标，而b则表示列下标，比如
A=[1,2,3;3,4,5] 则A([1,2],[2,3])=[2,3;4,5]
故最终结果有d(a)*d(b)个;d(a)表示a的维数。

模拟运动模型

clear all;
f=checkerboard(8);%生成一个8格的板子
g=pixeldup(f,8);

PSF = fspecial('motion',7,45);
gb=imfilter(f,PSF,'circular');
noise=imnoise(zeros(size(f)),'gaussian',0,0.001);
gbn=gb+noise;
myImshow(f,'拓展图');myImshow(gb,'运动后图像');myImshow(noise,'噪声图像');
myImshow(gbn,'运动后噪声图像');myImshow(0);

untitled.png

维纳滤波（最小均方误差误差）

用来做图像复原。其中H(u,v)是降指模型的点扩散函数

Paste_Image.png
Paste_Image.png