matlab spfilt,图像恢复

图像恢复

图像降指与恢复过程

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