chap4 图像复原 part2

最新推荐文章于 2022-03-08 14:47:29 发布

_helen_520

最新推荐文章于 2022-03-08 14:47:29 发布

阅读量568

点赞数

分类专栏：数字图像处理

本文链接：https://blog.csdn.net/haronchou/article/details/78071926

版权

数字图像处理专栏收录该内容

11 篇文章 3 订阅

订阅专栏

Reference：《数字图像处理 MATLAB版》冈萨雷斯

滤波时可用的所有工具函数：

g = imfilte(f, w, filtering_mode, boudary_options,size_options)要自己设计滤波器矩阵 w.
g = colfit(f,[m,n],'sliding',fun) 自己设计滤波器处理函数fun，同时注意A的元素排列规则

w = fspecial('type',parameters) 生成一些滤波器（laplace、均值、高斯、prewitt、LoG等）
冈萨雷斯设计的函数f = spfilt(g, type, varargin)

比较各个函数：imfilter也可以实现几何平均滤波等，比colfit要快得多，且内存需求也比colfit小的多。但若函数无法用imfilter实现的时候，colfilt还是非线性滤波最好的选择。
这里写图片描述

1. 非线性空间滤波

nlfilter 与 colfilt: nlfilt直接执行二维操作，而colfilt执行列操作，虽然占用内存，但速度加快。

几何平均 $\hat{f}(x,y) ={(\prod_{(s,t)\in{S_{xy}}}{g(s,t)})}^\frac{1}{mn}$

%以几何平均滤波为例 
gmean = @(A)prod(A,1)^(1/size(A,1));%其中prod(A),A为向量，则返回元素乘积；A为矩阵，则返回每列乘积组成的一个行向量；prod(A,dim)计算A中由dim指定方向的乘积。
f = padarray(f,[m,n],'replicate');%需要对图像填充；[m,n]表示填充多少行多少列；
g = colfit(f,[m,n],'sliding',@gmean);
%删除原来的填充
[M,N]=size(f);
g = g((1:M,(1:N));

fp = padarray(f,[r,c],method,direction);%method='symmetric' 'replicate' 'circular'; direction = 'pre' 'post' 'both'
%colfilt的使用说明
g = colfit(f,[m,n],'sliding',fun);%[m,n]指的是滤波器模板的大小；A=mn*MN，每列为MN每个像素的领域像素，fun对A的数据进行处理；如上面的几何平均函数gmean.同时将计算结果写成1*MN的行向量，再由colfilt函数将这个行向量变成M*N的图像大小；不过对于填充后的图像，要手动剪裁到原来大小

2. 工具箱里的标准空间滤波器

（1）线性空间滤波器 `w = fspecial('type',parameters)`

补充： LaTeX表示矩阵
高斯函数低通滤波器： $h(x,y)=e^{-\frac{x^2+y^2}{2\sigma^2}}$ ，对于3*3的滤波器： $w_1=h(-1,-1),w_2=h(-1,0)$ 等等

一阶微分算子： $\frac{\partial f}{\partial x}=f(x+1)-f(x)$ 对于两个元素的： $\nabla f=f(x+1，y)+f(x,y+1)-2f(x,y)$
二阶微分算子： $\frac{\partial^2 f}{\partial^2 x}=f(x+1)+f(x-1)-2f(x)$
laplace算子： $\nabla^2f(x,y)=f(x+1,y)+f(x-1,)+f(x,y+1)+f(x,y-1)-4f(x,y)$
laplace算子增加一个alpha: $\begin{matrix} \frac{\alpha}{1-\alpha}&\frac{1-\alpha}{1+\alpha}&\frac{\alpha}{1+\alpha} \\\frac{1-\alpha}{1+\alpha} & \frac{-4}{1+\alpha}& \frac{1-\alpha}{1+\alpha}\\ \frac{\alpha}{1+\alpha}&\frac{1-\alpha}{1+\alpha}&\frac{\alpha}{1+\alpha}\end{matrix}$ 可以对增强的结果进行精细的调整

一阶微分算子——梯度： $\nabla f =[g_x,g_y]'= [\frac{\partial f}{\partial x} \frac{\partial f}{\partial x}]'$ 则梯度的长度为： $M(x,y)=mag(\nabla f)=\sqrt{g_x^2+g_y^2}\approx |g_x|+|g_y|$ 由一阶微分算子可得: $g_x=z_8-z_5$ $g_y=z_6-z_5$ .
Robert提出微分交叉算子： $g_x=z_9-z_5, g_y=z_8-z_6$
Sobel算子： $M(x,y)\approx |g_x|+|g_y|=|(z_7+2z_8+z_9)-(z_1+2z_2+z_3)|+|(z_3+2z_6+z_9)-(z_1+2z_4+z_7)|$

例如：Robert算子 $g_x=\lgroup\begin{matrix} -1&0\\0&1\end{matrix}\rgroup$ $g_y=\lgroup\begin{matrix}0&-1\\1&0\end{matrix}\rgroup$
Sobel算子 $g_x=\begin{matrix}-1&-2&-1\\0&0&0\\1&2&1\end{matrix}$ $g_y=\begin{matrix}-1&0&-1\\-2&0&2\\-1&0&1\end{matrix}$

（2）fspecial 函数支持的空间滤波器如下：

这里写图片描述

(3) 非线性空间滤波器——统计排序滤波器

g = ordfilt2(f, order, domain) 使用领域的排序集合的第oder个元素去代替f中的每个元素。

g = ordfilt2(f, 1, ones(m,n)) 最小滤波器；g = ordfilt2(f, mn, ones(m,n)) 最大滤波器
g = ordfilt2(f, （m*n+1)/2, ones(m,n)) 中值滤波器
专门的中值滤波函数：g = medfilt2(f, [m,n], padopt) ;padopt为边界填充选项：’zeros’, ‘symmetric’, ‘indexed’ 表示若f是double类，用1填充；否则，用0填充。

median(-1:1)可用于计算这个区间的中值是多少

中值滤波器对于消除椒盐噪声的影响是很有效的！！！！！！！

3. 空间噪声滤波（chap4.3)

补充知识2 的表格中所列的滤波器用 spfilt函数来实现 split源代码见末尾

%% part2 非线性空间滤波——空间噪声滤波
clc;clear
f = imread('Fig0318(a)(ckt-board-orig).tif');
subplot(3,3,1);imshow(f)title('原始图像') 

[M,N] = size(f);
R = imnoise2('salt & pepper',M,N,0.1,0);
c = find(R == 0);
gp = f;
gp(c) = 0;
subplot(3,3,2);imshow(gp)title('被概率为0.1的胡椒噪声污染的图像')

R = imnoise2('salt & pepper',M,N,0,0.1);
c = find(R == 1);
gs = f;
gs(c) = 255;
subplot(3,3,3);imshow(gs)title('被概率为0.1的盐粒噪声污染的图像')

fp = spfilt(gp,'chmean',3,3,1.5);
subplot(3,3,4);imshow(fp)
title('用阶为Q=1.5的3*3反调和滤波器对[被概率为0.1的胡椒噪声污染的图像]滤波的结果')

fs = spfilt(gs,'chmean',3,3,-1.5);
subplot(3,3,5);imshow(fs)
title('用阶为Q=-1.5的3*3反调和滤波器对[被概率为0.1的盐粒噪声污染的图像]滤波的结果')

fpmax = spfilt(gp,'max',3,3);
subplot(3,3,6); imshow(fpmax)
title('用3*3最大滤波器对[被概率为0.1的胡椒噪声污染的图像]滤波的结果')

fsmin = spfilt(gs,'min',3,3);
subplot(3,3,7);imshow(fsmin)
title('用3*3最小滤波器对[被概率为0.1的盐粒噪声污染的图像]滤波的结果')

补充知识1：函数句柄——详见Reference

作用:一种数据类型，可作为参数调用。分为两类：一种命名函数句柄；一种匿名函数句柄。

%%命名句柄
f = @sin;
y = f(pi/4);%%等价于y = sin(pi/4)
%%匿名句柄
g = @(x)x.^2;
r = @(x,y)sqrt(x.^2 + y.^2);%不管是一个参数，还是两个参数，都可以

补充知识2：空域滤波器的表达式

这里写图片描述

补充3：split源代码

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: （M*N为滤波器模板大小）
%
%     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 修正的alpha平均滤波
%                                        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);