2017-03-26 12:18:38 u013288466 阅读数 4512
  • OpenCV3.2 Java图像处理视频学习教程

    OpenCV3.2 Java图像处理视频培训课程:基于OpenCV新版本3.2.0详细讲述Java OpenCV图像处理部分内容,包括Mat对象使用、图像读写、 基于常用核心API讲述基本原理、使用方法、参数、代码演示、图像处理思路与流程讲授。主要内容包括opencv像素操作、滤波、边缘提取、直线与圆检测、形态学操作与分水岭、图像金子塔融合重建、多尺度模板匹配、opencv人脸检测、OpenCV跟Tomcat使用实现服务器端图像处理服务。

    4137 人正在学习 去看看 贾志刚

眼底视网膜血管增强方法(一)匹配滤波

眼底是人体中唯一可以无创伤地观察到血管的地方,而血管的形态与人体疾病又有着重要的联系,因此眼底图像在疾病筛查中有着重要的意见。通过观察视网膜血管形状可以很方便地诊断出糖尿病、青光眼等病状。眼底血管的自动分割是计算机辅助诊断的一个重要手段,它可以极大的减少人力劳动,为广大患者提供方便快捷的医疗诊断
由于血管与背景的底对比度,使得进行血管分割的难度增大,因此,血管增强是血管分割中常用的预处理方法。视网膜血管增强有很多方法,像形态学增强、边缘检测增强、滤波增强等,其中最为经典的当为CHAUDHURI1的匹配滤波方法。

信号处理中的匹配滤波

匹配滤波来源于通信系统中的一维信号处理。在通信系统中,匹配滤波接收器能够使得接收信号的输出信噪比最大。设滤波器的传递函数为H(f),冲击响应为h(t),滤波器输入码元s(t),的持续时间为Ts,则输入信号r(t)

r(t)=s(t)+n(t)0tTs

其中n(t)为高斯白噪声。
则根据Schwarz定理得,当
H(f)=kS(f)ej2πft0
时,得到最大的输出信噪比2E/n,即匹配时要求h(t)=s(t).

眼底图像的匹配滤波

匹配原理

眼底图像的血管成管状,有很好的形态学特性。血管的横跨面的灰度分布成高斯形状,而背景的灰度又基本一至,如图一所示。因此,我们可以把血管的横截面(一维,即垂直于血管的一段线段)看着是掺杂着高斯白噪声(背景点)的一维信号(血管点)。根据信号处理的匹配滤波原理,我们可以选取一传递函数与血管分布一致的(高斯型)滤波器来得到最大信噪比的输出,即可增强血管。
图一 血管中心灰度分布图
图一 血管中心灰度分布图

匹配滤波器的构造

要使得滤波器与血管相匹配,则滤波器有两个重要的参数需要考虑,一是尺度,二是方向。

尺度的匹配

从图一可以看出,不同大小的血管的分布的尺度(对应高斯函数的带宽σ)是不一样的。为了达到更好的匹配效果,同时能增强大血管和小血管,文章选取的4个不同尺度(2,2,22,4)的滤波器来进行滤波。(能不能构造一种尺度也自适应的匹配滤波?)

方向的匹配

由于滤波器的方向必须垂直于血管的方向,而眼底中血管的方向又是随机伸展的。因此,文章构造的12个方向的滤波器,在每个像素点都进行12次不同方向的滤波,选取最大响应的一个作为最终的响应输出。(能不能构造一种方向也自适应的匹配滤波?)同时,血管是现状的,在一定的长度范围内,他们的方向是相同的,因为我们可以对一yl2的小块进行同时滤波,这样可以提高效率。

滤波器的表达式

根据上述,我们可以得到滤波器的数学表达式

K(x,y)=exp(x22δ2)yL2

实验结果

原图像 原图像

原程序

    sigma=2;
    yLength=9;
    direction_number=12;
    MF = MatchFilter(img, sigma, yLength,direction_number);
    MF(mask==0) = 0;
    MF = normalize(double(MF));
    % Adding to features
    features = MF;

程序包下载



  1. SUBHASIS CHAUDHURI, STUDENT MEMBER,”Detection of Blood Vessels in Retinal Images Using
2017-10-12 13:01:27 piaoxuezhong 阅读数 8071
  • OpenCV3.2 Java图像处理视频学习教程

    OpenCV3.2 Java图像处理视频培训课程:基于OpenCV新版本3.2.0详细讲述Java OpenCV图像处理部分内容,包括Mat对象使用、图像读写、 基于常用核心API讲述基本原理、使用方法、参数、代码演示、图像处理思路与流程讲授。主要内容包括opencv像素操作、滤波、边缘提取、直线与圆检测、形态学操作与分水岭、图像金子塔融合重建、多尺度模板匹配、opencv人脸检测、OpenCV跟Tomcat使用实现服务器端图像处理服务。

    4137 人正在学习 去看看 贾志刚

       视网膜血管分割是眼科计算机辅助诊断和大规模疾病筛查系统的基础,当眼器官发生视觉疾病的时候,网膜血管的直径、颜色和弯曲程度等就会出现异常,科医生可以据此作出诊断。医疗上对血管图像进行分割常用的方法有:基于血管跟踪的方法、基于匹配滤波的方法、基于形态学处理的方法、基于形变模型的方法和基于机器学习的方法等 。本篇将介绍基于匹配滤波算法。

匹配滤波器(matched filter)算法

原理部分:

 

      匹配滤波算法最早见于论文【1】《Detection of  Blood Vessels in Retinal Images Using Two-Dimensional Matched Filters

 

》中,后来又有许多针对该算法的不足进行改进的算法出现。

      匹配滤波器是指经过滤波后,使得滤波器输出端的信号瞬时功率与噪声平均功率的比值(即是信噪比SNR)最大,当有用信号与噪声同时进入滤波器时,有用信号在某一瞬间出现尖峰值,而噪声信号受到抑制。

      论文【1】对眼底视网膜上的血管特点进行分析,基于这些特点设计匹配滤波器,主要参考特点:血管的宽度在较小的范围内变动,血管壁的两条线是近似平行的,血管有方向,血管横切线上的灰度曲线如下,近似满足高斯分布:

        要理解论文中的公式推导,则需要对匹配滤波有更深层次的理解和掌握。首先,匹配滤波器也是一种在某一准则条件下的最优设计,其对应的最优的准则是输出信噪比(SNR)最大,而且前提条件是噪声信号为白噪声。论文得到的匹配滤波器的表达式为:H(f)=S*(f),即匹配滤波器的频率响应是输入信号频率响应的共轭。从物理直观上理解:

(1)从幅频特性来看,匹配滤波器和输入信号的幅频特性一致。在信号越强的频率点,滤波器的放大倍数也越大;在信号越弱的频率点,滤波器的放大倍数越小。这就是信号处理中的“马太效应”,匹配滤波器是让有用信号尽可能通过。白噪声的功率谱在各个频率点都一样,这种情况下,有用信号尽可能通过,噪声则相反。
(2)从相频特性上看,匹配滤波器的相频特性和输入信号则相反。通过匹配滤波器后,信号的相位为0,正好能实现信号时域上的相干叠加。而噪声的相位是随机的,只能实现非相干叠加。这样在时域上保证了输出信噪比的最大。
       在信号的幅频特性与相频特性中,幅频特性表征了频率特性,频特性表征了时间特性。匹配滤波器无论是从时域还是从频域,都充分保证了信号尽可能大的通过,噪声尽可能小的通过。

       回到论文【1】中的描述,将血管想象成一小段一小段的平行区域的组合,设定长度为L,宽度为3σ,然后用一个高斯曲线来模拟血管横切面上的灰度曲线, 匹配滤波器如下:

原文中讲到:“When the concept of matched filter is extended to two- dimensional images, it must be appreciated that a vessel may be oriented at any angle θ (0 <=θ<=pi ). The matched filter s ( t ) will have its peak response only when it is aligned at an angle θ+- pi / 2 . ” (不是很理解为什么,有知道的请告知?)

由于眼底中血管是有方向的,并且血管的方向是随机伸展的,而滤波器在垂直于血管的方向上出现峰值。所以设计滤波器从0度到180度每15度旋转一次,得到12个方向上的滤波器,然后分别进行卷积,选取最大响应的一个作为最终的响应输出。

旋转矩阵为:

算法实现中用到的核函数及掩码设置参见原文,怕翻的不精准(捂脸~)

算法实现:

参考1中的代码试过效果不是很好,可能是参数需要调试,这里附录下:

 

%% matched filter2,测试程序
clc,close all,clear all;
img=imread('../image/1.jpg');
if length(size(img))==3
    img=rgb2gray(img);
end
[g,bg]=matchedFilter2(img);

 

 

function [g,bg]=matchedFilter2(f)
 
f=double(f);
subplot(1,2,1);
imshow(uint8(f));
title('origin image');
% mean filter
f=medfilt2(f,[5,5]);
% f=medfilt2(f,[21 1]);
% f=medfilt2(f,[1,7]);
% 参数
os=12;  % 角度的个数
sigma=2;
tim=3;
L=9;
t=180; % 全局阈值,需要多次尝试
 
thetas=0:(os-1);
thetas=thetas.*(180/os);
N1=-tim*sigma:tim*sigma;
N1=-exp(-(N1.^2)/(2*sigma*sigma));
N=repmat(N1,[2*floor(L/2)+1,1]);
r2=floor(L/2);
c2=floor(tim*sigma);
[m,n]=size(f);
RNs=cell(1,os);  % rotated kernals
MFRs=cell(1,os); % filtered images
g1=f;
 
% matched filter
for i=1:os
    theta=thetas(i);
    RN=imrotate(N,theta);
    %去掉多余的0行和零列
    RN=RN(:,any(RN));
    RN=RN(any(RN'),:);
    meanN=mean2(RN);
    RN=RN-meanN;
    RNs{1,i}=RN;
    MFRs{1,i}=imfilter(f,RN,'conv','symmetric');
end
 
% get the max response
g=MFRs{1,1};
for j=2:os
    g=max(g,MFRs{1,j});
end
bg=g<t;

subplot(1,2,2);
imshow(bg);
title('filterd image');
end

上述代码中的三个参数sigma,L,T的选择需要针对不同的应用去尝试,如果是视网膜血管图像,参数可以设置为2,9,3。全局阈值t的选择,可以多试几次找最好的。测试结果如下:

 

我重点测试了参考3的算法实现,测试代码如下,整个工程文件见后面链接:

 

 

%% 测试函数,匹配滤波算法血管分割
clc,close all,clear all;
sigma=2;
yLength=10;
direction_number=12;
img=imread('../image/1.bmp');
if length(size(img))==3
    img=rgb2gray(img);
end
MF = MatchFilter(img, sigma, yLength,direction_number);
mask=[0 0 0 0 0;
            0 1 1 1 0;
            0 1 1 1 0;
            0 1 1 1 0;
            0 0 0 0 0;];
MF(mask==0) = 0;
MF = normalize(double(MF));
% Adding to features
features = MF;
subplot(1,2,1);
imshow(img);title('origin image');
subplot(1,2,2);
imshow(features);title('filted image');

效果图如下:

 

全部工程文件下载链接

算法优缺点:

滤波之后图像中的血管结构得到增强,但是该方法可能会丢失血管分叉点和细小的血管;

参考:

 

  • http://www.cnblogs.com/naniJser/archive/2012/12/09/2810551.html
  • http://blog.csdn.net/chiooo/article/details/43564069
  • http://blog.csdn.net/u013288466/article/details/66473284

 

2012-02-27 13:28:13 lj695242104 阅读数 25772
  • OpenCV3.2 Java图像处理视频学习教程

    OpenCV3.2 Java图像处理视频培训课程:基于OpenCV新版本3.2.0详细讲述Java OpenCV图像处理部分内容,包括Mat对象使用、图像读写、 基于常用核心API讲述基本原理、使用方法、参数、代码演示、图像处理思路与流程讲授。主要内容包括opencv像素操作、滤波、边缘提取、直线与圆检测、形态学操作与分水岭、图像金子塔融合重建、多尺度模板匹配、opencv人脸检测、OpenCV跟Tomcat使用实现服务器端图像处理服务。

    4137 人正在学习 去看看 贾志刚

描述:均值滤波器是图像处理中一种常见的滤波器,它主要应用于平滑噪声。它的原理主要是利用某像素点周边像素的平均值来打到平滑噪声的效果。

常用的均值核如下图所示:


            


图像滤波器操作实际上就是模板操作,对于模板操作我们应该注意边界问题:

什么是边界问题?

对于边界问题就是当图像处理边界像素的时候,卷积核与图像使用区域不能匹配,计算出现问题。


处理方法:

1、忽略边界像素,即丢掉不能匹配的像素

2、保留边界像素,即复制源图像的不能匹配的边界像素到输出图像



Code:

  /**
   * Calculates the mean of a 3x3 pixel neighbourhood (including centre pixel).
   *
   * @param input the input image 2D array
   * @param kernel the kernel 2D array
   * @param w the image width
   * @param h the image height
   * @param x the x coordinate of the centre pixel of the array
   * @param y the y coordinate of the centre pixel of the array
   * @return the mean of the 9 pixels
   */ 
  public static int meanNeighbour(int [][] input, int [][] kernel,
			  int w, int h, int x, int y) {

    int sum = 0;
    int number = 0;
    for(int j=0;j<3;++j){
      for(int i=0;i<3;++i){
	if((kernel[i][j]==1) && 
	   ((x-1+i)>=0) && ((y-1+j)>=0) && ((x-1+i)<w) && ((y-1+j)<h) ){
	  sum = sum + input[x-1+i][y-1+j];
	  ++number;
	}
      }
    }
    if(number==0) return 0;
    return (sum/number);
  }


  /**
   * Takes an image in 2D array form and smoothes it according to the kernel.
   * @param input the input image
   * @kernel the kernel 1D array
   * @param width of the input image
   * @param height of the output image
   * @param iterations to be performed
   * @return the new smoothed image 2D array
   */
  public static int [][] smooth(int [][] input, int [][] kernel,
				int width, int height, int iterations){
    int [][] temporary = new int [width][height];
    int [][] outputArrays = new int [width][height];
    temporary = (int [][]) input.clone();
    for (int its=0;its<iterations;++its){
      for(int j=0;j<height;++j){
	for(int i=0;i<width;++i){
	  outputArrays[i][j] = meanNeighbour(temporary,kernel,
					     width,height,i,j);
	}
      }
      for(int j=0;j<height;++j){
	for(int i=0;i<width;++i){
	  temporary[i][j]=outputArrays[i][j];
	}
      }
    }
    return outputArrays;
  }


Input Image:




Output Image:




总结:均值滤波器就是为了平滑周边的效果,不会因为图像上突出的点而感到难受,达到一种“Softened”的效果。

2019-06-27 14:08:54 weixin_41683113 阅读数 124
  • OpenCV3.2 Java图像处理视频学习教程

    OpenCV3.2 Java图像处理视频培训课程:基于OpenCV新版本3.2.0详细讲述Java OpenCV图像处理部分内容,包括Mat对象使用、图像读写、 基于常用核心API讲述基本原理、使用方法、参数、代码演示、图像处理思路与流程讲授。主要内容包括opencv像素操作、滤波、边缘提取、直线与圆检测、形态学操作与分水岭、图像金子塔融合重建、多尺度模板匹配、opencv人脸检测、OpenCV跟Tomcat使用实现服务器端图像处理服务。

    4137 人正在学习 去看看 贾志刚

术语:空间域

空间域指图像平面本身,这是相对于变换域而言的,变换域的图像处理首先把空间域变换到变换域,在变换域中处理图像,再通过反变换,把处理结果返回到空间域。

空间域的图像处理直接对图像中的像素进行操作,主要分为灰度变换和空间滤波。

空间滤波也称领域处理或空间卷积。

空间域处理的表示:

g(x,y)=T[f(x,y)] g(x,y)=T[f(x,y)]

  • f(x,y)是输入图像,g(x,y) 是输出图像,T 是在点 (x,y) 的邻域上定义的关于 f 的一种算子。

  • T 可作用于单幅图像,也可作用于图像集合,如为了降低噪声而对图像序列执行逐像素求和操作。

  • 领域通常取为以 (x,y) 为中心的矩形领域,大小为 (2a+1)×(2b+1)(2a+1)\times(2b+1)

  • 空间滤波的过程:领域中心从图像原点开始,从左到右,从上到下移动,输出图像在每一个空间坐标 (x,y)
    处的像素值就是 T 作用在 f 的以 (x,y) 为中心的领域上所得的值。

    note:当邻域中心充分靠近图像边界时,部分领域会位于图像外部,此时需要做一些技术处理,主要有以下两种方法:

    • 忽略外侧邻点
    • 用0或者其他指定灰度值填充图像边缘
      3.1

术语:空间滤波器

领域与预定义的操作 T 一起称为空间滤波器,也称空间掩模,核,模板或窗口。

灰度变换(点处理技术)

当领域取为最小领域 1×11\times 1 时,g 仅取决于 (x,y) 处的 f 值与算子 T .此时空间域处理可以简单得表示成:
s=T(r) s=T(r)

  • r,s 分别表示 f,g 在任意点 (x,y) 的灰度

灰度变换函数

线性函数

  • 恒等变换
    不改变图像,仅出于理论完整性考虑

  • 反转变换

    • 表达式:
      s=L1r s=L-1-r

    • 此处假设灰度级范围是[0 , L-1]

    • 作用:可以得到等效的照片底片,适用于增强嵌入在暗区域中白色或灰色细节,特别是黑色占主导地位时。

    • 代码实现:

       function g = imreverse(f)
       %IMREVERSE computes the negative of an image f
       
       % Convert the image to uint8_type
       f_uint8 = im2uint8(f);
       g = 255 - f_uint8;
      
    • 效果展示:
      input:

       >> f = imread('Fig0304(a)(breast_digital_Xray).tif');
       >> g = imreverse(f);
       >> figure
       >> subplot(1,2,1);imshow(f),title('原图');
       >> subplot(1,2,2);imshow(g),title('图像反转');
      

      3.12

对数函数

  • 通用形式:
    s=clog(1+r) s = clog(1+r)
  • c是常数,假设 r0r\ge 0
  • 图像:
    3.18
  • 从对数曲线的形状可以看出,该变换将较窄范围的暗色输入值映射为较宽范围的输出值,将较宽范围的亮色输入值映射为较窄范围的输出值。
  • 代码实现:
     function g = imlog(f)
     %IMLOG takes the logarithm transform of an image f
     
     f1 = mat2gray(im2double(f));
     g = log2(1+f1);
    
  • 效果展示:
     >> f = imread('Fig0316(4)(bottom_left).tif');
     >> g = imlog(f);
     >> figure
     >> subplot(1,2,1);imshow(f),title('原图');
     >> subplot(1,2,2);imshow(g),title('log');
    
    3.25
  • 作用:
    • 放大低灰度值区域的细节
    • 压缩动态范围
      例:当因为数值范围过大,导致图像显示效果不佳时,可以利用对数运算将大范围压缩到小范围这一性质,做对数变换后,再显示图像。
      以下原图是范围在0到1.5\times 10^6 的傅里叶频谱,频谱中的低值损失严重。做对数变换后,效果好很多。
       f = imread('Fig0305(a)(DFT_no_log).tif');
       g = im2uint8(mat2gray(log(1+double(f))));
       figure
       subplot(1,2,1);imshow(f),title('原图');
       subplot(1,2,2);imshow(g),title('log');
      
      3.19

伽马变换

  • 基本形式:
    s=crγ s = cr^{\gamma}

    • c,γ\gamma 是正常数
    • γ&gt;1\gamma &gt; 1 时,该变换将较窄范围的暗色输入值映射为较宽范围的输出值,将较宽范围的亮色输入值映射为较窄范围的输出值。
      γ&lt;1\gamma&lt;1 时,效果相反
    • 不同 γ\gamma 值下,曲线的形状:
      3.13
  • 代码实现:
    M文件:

     function g = imgamma(f, gamma)
     %IMGAMMA takes a gamma transform of an image f
     
     % First, we normalizes the pixel value of the image f to [0,1]
     f1 = double(f);
     f2 = mat2gray(f1);
     g = f2.^gamma;
    

    例:
    input:

     >> f = imread('Fig0308(a)(fractured_spine).tif');
     >> g = imgamma(f,0.3);
     >> figure
     >> subplot(1,2,1);imshow(f),title('原图');
     >> subplot(1,2,2);imshow(g),title('gamma = 0.3');
    

    output:
    3.15

  • 应用:

    • 伽马校正:用于图像获取,打印和显示的各种设备自身根据伽马变换产生响应。显示系统的这一特性会使显示的图像和原图像相比偏暗或偏亮。这时根据设备的伽马值做相应的伽马变换,称为伽马校正,使得设备能正确显示原图。
      例:下图是监视器上显示图像的伽马校正过程。因为CRT有一个灰度-电压响应,相当于显示图像时自动做了一个 γ\gamma 值在1.8到2.5的伽马变换,使得监视器显示的图像偏暗。在显示图像前,先对原图做一个伽马值等于 1γ\frac{1}{\gamma} 的伽马变换作为校正,这样监视器就能正确显示原图像了。
      3.14

    • 对比度增强:
      例:有时图像偏暗,可以使用 γ&lt;1\gamma &lt;1 的伽马变换。为了使得显示效果最佳,通常要尝试不同的 γ\gamma 值。
      γ\gamma 越接近1,图像越接近原图;γ\gamma 越接近0,原图暗处的细节被增强,但对比度降低,图像容易被"冲淡"。

       f = imread('Fig0308(a)(fractured_spine).tif');
       gamma = [0.6 0.4 0.3];
       figure
       subplot(1,4,1);imshow(f),title('原图');
       for i = 1:3
           g = imgamma(f,gamma(i));
           subplot(1,4,i+1);imshow(g),title(gamma(i));
       end
      

      3.16
      γ=0.3\gamma = 0.3 时,黑色背景已出现很多灰点,γ=0.4\gamma=0.4 效果最好

      例:有时图像偏亮,可以使用 γ&gt;1\gamma &gt;1 的伽马变换。为了使得显示效果最佳,通常要尝试不同的 γ\gamma 值。
      γ\gamma 越接近1,图像越接近原图,但原图偏亮,图像有"冲淡"的外观;γ\gamma 越接近无穷,原图越暗,太暗会丢失细节。

       f = imread('Fig0309(a)(washed_out_aerial_image).tif');
       gamma = [3 4 5];
       figure
       subplot(1,4,1);imshow(f),title('原图');
       for i = 1:3
           g = imgamma(f,gamma(i));
           subplot(1,4,i+1);imshow(g),title(gamma(i));
       end
      

      3.17
      γ=5\gamma = 5 时,图像已经偏暗,暗处细节已经丢失,γ=4\gamma=4 效果最好

对比度拉伸变换函数

  • 基本形式:
    s=T(r)=11+(m/r)E s=T(r)=\frac{1}{1+(m/r)^E}

    • 参数 E,m 控制函数的形状
  • 图像:
    3.24
    可以看出:这个变换将小于 m 的输入灰度级压缩到靠近0的窄区域,将大于 m 的输入灰度级压缩到靠近1的窄区域,使得输出具有高对比度。E越小,压缩效果越小;E越大,压缩能力越大,图像越向二值图像靠近。

  • 代码实现

     function g = imContrasStretch(f, m, E)
     %IMCONTRASSTRETCH takes the ContrasStretch tansform of an image f
     % m is the x-value of the inflection point
     % E determines the steepness of the curve
     
     f1 = mat2gray(im2double(f));
     g = 1./(1+(m./(f1+eps)).^E);
    
  • 效果展示:

     >> f = imread('rose_512.tif');
     >> g = imContrasStretch(f,0.5,20);
     >> figure
     >> subplot(1,2,1);imshow(f),title('原图');
     >> subplot(1,2,2);imshow(g),title('对比度拉伸');
    

    3.26

阈值处理函数

  • 基本形式:
    s={0r&lt;m1rm s=\left\{ \begin{array}{rcl} 0 &amp; &amp; {r&lt;m}\\ 1 &amp; &amp; {r\ge m}\\ \end{array} \right.
    • 将输入图像变为二值图像,像素值小于 m 变换为0,大于或等于 m 变换为1
  • 图像:
    3.23
  • 代码实现:
     function g = myIm2bw(f,T)
     %MYIM2BW turns an image f into a logical image
     % T is the threshold that belongs to [0,1]
     
     f1 = mat2gray(im2double(f));
     g = (f1 >= T);
    
  • 效果展示及与MATLAB内置函数im2bw对比
     >> f = imread('rose_512.tif');
     >> g1 = myIm2bw(f,0.5);
     >> g2 = im2bw(f,0.5);
     >> figure
     >> subplot(1,3,1);imshow(f),title('原图');
     >> subplot(1,3,2);imshow(g1),title('myIm2bw');
     >> subplot(1,3,3);imshow(g2),title('Im2bw');
    
    3.27
  • 应用:图像分割

分段线性变换

  • 优点:可以任意复杂,满足各种需求。
  • 缺点:要求用户寻找合适的函数
  • 应用:
    • 对比度拉伸:扩展图像灰度级动态范围

      • 例:构造以下灰度变换(像素值归一化为[0,1])
        s={(s1/r1)rr&lt;r1s2s1r2r1(rr1)+s1r1rr21s21r2(rr2)+s2r&gt;r2 s=\left\{ \begin{array}{rcl} (s_1/r_1)\cdot r&amp; &amp; {r&lt;r_1}\\ \frac{s_2-s_1}{r_2-r_1}(r-r_1)+s_1 &amp; &amp; {r_1\le r\le r_2}\\ \frac{1-s_2}{1-r_2}(r-r_2)+s_2&amp; &amp; {r&gt;r_2} \end{array} \right.
    • 图像:
      3.41

    • 代码实现

       function g = imcontras(f,r1,s1,r2,s2)
       %IMCONTRAS take the piecewise linear transform of an image f
       % r1<r2,s1<s2,(r1,s1) and (r2,s2) are turning points
       
       ff = mat2gray(im2double(f));
       [M,N] = size(f);
       g = zeros(M,N);
       for i = 1:M
           for j = 1:N
               if ff(i,j)<r1
                   g(i,j) = (s1/r1)*ff(i,j);
               elseif ff(i,j)<r2
                   g(i,j) = ((s2-s1)/(r2-r1))*(ff(i,j)-r1)+s1;
               else
                   g(i,j) = ((s2-1)/(r2-1))*(ff(i,j)-r2)+s2;
               end
           end
       end
      
    • 效果展示

        >> f = imread('Fig0320(2)(2nd_from_top).tif');
        >> g = imcontras(f,0.3,0.125,0.6,0.875);
        >> figure
        >> subplot(1,2,1);imshow(f),title('原图');
        >> subplot(1,2,2);imshow(g),title('分段线性变换后');
      

      3.42

    • 灰度级分层
      作用:突出特定灰度范围
      基本方法1:将感兴趣的灰度值变为1或者接近1,其他值变为0或者接近0

      • 图像:
        3.43
      • 代码实现:
         function g = imlayered(f,r1,s1,r2,s2)
         %IMCONTRAS take the piecewise linear transform of an image f
         % r1<r2,(r1,s1) and (r2,s2) are turning points
         % s1 is close to 0, s2 is close to 1
         
         ff = mat2gray(im2double(f));
         [M,N] = size(f);
         g = zeros(M,N);
         for i = 1:M
             for j = 1:N
                 if ff(i,j)<r1 || ff(i,j)>r2
                     g(i,j) = s1;
                 else
                     g(i,j) = s2;
                 end
             end
         end
        
      • 效果展示
         f = imread('Fig0235(c)(kidney_original).tif');
         g = imlayered(f,0.7,0,0.9,1);
         figure
         subplot(1,2,1);imshow(f),title('原图');
         subplot(1,2,2);imshow(g),title('灰度分层后');
        
        3.44
        基本方法2:将感兴趣的灰度值变为1或者接近1,其他值不变
      • 图像:
        3.45
      • 代码实现
         function g = imlayered2(f,r1,r2,s)
         %IMCONTRAS take the piecewise linear transform of an image f
         % r1<r2
         % s is close to 1
         ff = mat2gray(im2double(f));
         [M,N] = size(f);
         g = zeros(M,N);
         for i = 1:M
             for j = 1:N
                 if ff(i,j)<r1 || ff(i,j)>r2
                     g(i,j) = ff(i,j);
                 else
                     g(i,j) = s;
                 end
             end
         end
        
      • 效果展示
         f = imread('Fig0235(c)(kidney_original).tif');
         g = imlayered2(f,0.6,0.9,1);
         figure
         subplot(1,2,1);imshow(f),title('原图');
         subplot(1,2,2);imshow(g),title('灰度分层后');
        
        3.46
    • 比特平面分层

      • 作用:突出特定比特

      术语:比特平面:对于一个8比特的图像,共8层比特平面,第 i(i=2,3,...,8)i(i=2,3,...,8) 层比特平面指像素值在 [2i1,2i1][2^{i-1},2^{i}-1], 第一层是 {0,1}\{0,1\}

      • 总结:低比特平面像素值范围小,包含的图像信息少,它们主要贡献更精细的细节
        高比特平面像素值范围大,包含了图像的主要信息
      • 代码实现
         function imBitLayer(f)
         %IMBITLAYER receive an image f and return its bitplanes
         
         ff = im2uint8(f);
         [M,N] = size(f);
         for k = 2:8
             g(:,:,k) = zeros(M,N);
             for i = 1:M
                 for j = 1:N
                     if ff(i,j)>=2^(k-1) && ff(i,j)<=2^k-1
                         g(i,j,k)=1;
                     end
                 end
             end
         end
         for i = 1:M
             for j = 1:N
                 if ff(i,j)==0 || ff(i,j)==1
                         g(i,j,1)=1;
                 end
             end
         end
         figure
         subplot(3,3,1);imshow(f),title('原图');
         for i = 2:9
             subplot(3,3,i);imshow(g(:,:,i-1)),title(i-1);
         end
        
      • 效果展示
         >> f = imread('Fig0314(a)(100-dollars).tif');
         >> imBitLayer(f)
        
        3.47

IPT提供的API

imadjust

  • 语法:
    g=imadjust(f,[low_in,high_in],[low_out,high_out],gamma) g=imadjust(f,[low\_in,high\_in],[low\_out,high\_out],gamma)

  • 输入图像应属于uint8,uint16或double类。输出图像与输入图像同类型。

  • [low_in,high_in] , [low_out , high_out] 都是[0,1]的子区间,若设置为[],则表示[0,1]
    它们与gamma共同决定了一个分段幂律变换,公式如下:
    不妨记:[ low_in , high_in] = [a,b] , [low_out , high_out]=[c,d] , gamma= γ\gamma
    g={cf&lt;a(dc)(faba)γ+cafbdf&gt;b g=\left\{ \begin{array}{rcl} c &amp; &amp; {f&lt;a}\\ (d-c)(\frac{f-a}{b-a})^{\gamma}+c &amp; &amp; {a\le f\le b}\\ d &amp; &amp; {f&gt;b} \end{array} \right.

  • 参数gamma指明了曲线形状。gamma<1, 映射加权至较高输出值 ;gamma>1,映射加权至较低输出值。gamma默认为1
    3.3

  • imadjust先将 f 归一化,再代入上述公式计算g ,最后根据 f 的类型将g的值还原为相应的范围,过程中遇到需要取整时就取整。

  • 特别地,high_out可以小于low_out,此时输出亮度将反转

  • imadjust可以实现图像反转,伽马变换和某些分段线性变换

  • 例:以下图片是一幅数字乳房X射线图像,显示出一处病灶,原图偏亮的病灶嵌在一大片黑色区域中,不利于分析病情,使用imadjust做图像增强处理。

       >> f = imread('Fig0304(a)(breast_digital_Xray).tif');
       >> g1 = imadjust(f,[0,1],[1,0]);
       >> g2 = imadjust(f,[0.5 0.75],[0 1]);
       >> g3 = imadjust(f,[],[],2);
       >> figure
       >> subplot(1,4,1);imshow(f),title('原始图像');
       >> subplot(1,4,2);imshow(g1),title('g1');
       >> subplot(1,4,3);imshow(g2),title('g2');
       >> subplot(1,4,4);imshow(g3),title('g3');
    

    输出:
    3.6

    • g1是原图的图像反转
    • g2是将原图位于[0.5 , 0.75]之间的亮度扩展到[0 , 1],用于强调图像中感兴趣的亮度区域。
    • g3是原图的伽马变换,压缩了低端,扩展了高端

note: 函数imcomplement也可以得到图像反转。

   >> f = imread('Fig0304(a)(breast_digital_Xray).tif');
   >> g = imcomplement(f);
   >> figure
   >> subplot(1,2,1);imshow(f),title('原图');
   >> subplot(1,2,2);imshow(g),title('负片');

输出:

3.8

stretchlim

用于自动设置imadjust的[ low_in , high_in]

  • 语法:
    Low_High=stretchlim(f,tol) Low\_High = stretchlim(f , tol)

    • 如果tol是二维向量[low_frac , high_frac],指定了图像低和高像素值饱和度的百分比。
    • 如果tol是标量,那么[low_frac , high_frac] = [tol, 1-tol]
    • tol默认为[0.01,0.99]
    • tol = 0,low_high = [min(f( : )) , max(f( : ))]
  • 例:

     >> f = imread('Fig0304(a)(breast_digital_Xray).tif');
     >> g = imadjust(f,stretchlim(f),[]);
     >> figure
     >> subplot(1,2,1);imshow(f),title('原图');
     >> subplot(1,2,2);imshow(g),title('饱和度水平2%');
    

    输出:

    3.7

灰度变换的一些实用M函数

处理可变数目的输入或输出

  • nargin: 返回输入到M函数的输入变量数
  • nargout: 返回M函数的输出变量数
  • nargchk:
    • 语法:
      msg=nargchk(low,high,number) msg=nargchk(low,high,number)
      • 检查传递的参数数量是否正确
      • number<low,返回“Not enough input parameters”;number>high, 返回“Too many input parameters”,number介于两者之间,返回空矩阵
  • varargin/varargout: 用于可变数量的参数,以单元数组的形式存储

intrans函数的实现

  • 代码:

      function g = intrans(f, varargin)
      %INTRANS Performs intensity (gray-level) transformations.
      %   G = INTRANS(F, 'neg') computes the negative of input image F.
      % 
      %   G = INTRANS(F, 'log', C, CLASS) computes C*log(1 + F) and
      %   multiplies the result by (positive) constant C. If the last two
      %   parameters are omitted, C defaults to 1. Because the log is used
      %   frequently to display Fourier spectra, parameter CLASS offers the
      %   option to specify the class of the output as 'uint8' or
      %   'uint16'. If parameter CLASS is omitted, the output is of the 
      %   same class as the input. 
      % 
      %   G = INTRANS(F, 'gamma', GAM) performs a gamma transformation on
      %   the input image using parameter GAM (a required input).  
      %
      %   G = INTRANS(F, 'stretch', M, E) computes a contrast-stretching
      %   transformation using the expression 1./(1 + (M./(F +
      %   eps)).^E).  Parameter M must be in the range [0, 1].  The default
      %   value for M is mean2(im2double(F)), and the default value for E
      %   is 4.
      %
      %   For the 'neg', 'gamma', and 'stretch' transformations, double
      %   input images whose maximum value is greater than 1 are scaled
      %   first using MAT2GRAY.  Other images are converted to double first
      %   using IM2DOUBLE.  For the 'log' transformation, double images are
      %   transformed without being scaled; other images are converted to
      %   double first using IM2DOUBLE.
      %
      %   The output is of the same class as the input, except if a
      %   different class is specified for the 'log' option.
      
      % Verify the correct number of inputs.
      error(nargchk(2, 4, nargin))
      
      % Store the class of the input for use later.
      classin = class(f);
      
      % If the input is of class double, and it is outside the range
      % [0, 1], and the specified transformation is not 'log', convert the
      % input to the range [0, 1].
      if strcmp(class(f), 'double') & max(f(:)) > 1 & ...
            ~strcmp(varargin{1}, 'log')
         f = mat2gray(f);
      else % Convert to double, regardless of class(f).
         f = im2double(f);
      end
      
      % Determine the type of transformation specified.
      method = varargin{1};
      
      % Perform the intensity transformation specified.    
      switch method
      case 'neg' 
         g = imcomplement(f); 
      
      case 'log'
         if length(varargin) == 1  
            c = 1;
         elseif length(varargin) == 2  
            c = varargin{2}; 
         elseif length(varargin) == 3 
            c = varargin{2}; 
            classin = varargin{3};
         else 
            error('Incorrect number of inputs for the log option.')
         end
         g = c*(log(1 + double(f)));
      
      case 'gamma'
         if length(varargin) < 2
            error('Not enough inputs for the gamma option.')
         end
         gam = varargin{2}; 
         g = imadjust(f, [ ], [ ], gam);
         
      case 'stretch'
         if length(varargin) == 1
            % Use defaults.
            m = mean2(f);  
            E = 4.0;           
         elseif length(varargin) == 3
            m = varargin{2};  
            E = varargin{3};
         else error('Incorrect number of inputs for the stretch option.')
         end
         g = 1./(1 + (m./(f + eps)).^E);
      otherwise
         error('Unknown enhancement method.')
      end
      
      % Convert to the class of the input image.
      g = changeclass(classin, g);
    
  • 例:

     f = imread('Fig0306(a)(bone-scan-GE).tif');
     g = intrans(f,'stretch',mean2(im2double(f)),0.9);
     figure
     subplot(1,2,1);imshow(f),title('原图');
     subplot(1,2,2);imshow(g),title('stretch');
    

    3.48

亮度标度的M函数:gscale

  • 代码:

     function g = gscale(f, varargin)
     %GSCALE Scales the intensity of the input image.
     %   G = GSCALE(F, 'full8') scales the intensities of F to the full
     %   8-bit intensity range [0, 255].  This is the default if there is
     %   only one input argument.
     %
     %   G = GSCALE(F, 'full16') scales the intensities of F to the full
     %   16-bit intensity range [0, 65535]. 
     %
     %   G = GSCALE(F, 'minmax', LOW, HIGH) scales the intensities of F to
     %   the range [LOW, HIGH]. These values must be provided, and they
     %   must be in the range [0, 1], independently of the class of the
     %   input. GSCALE performs any necessary scaling. If the input is of
     %   class double, and its values are not in the range [0, 1], then
     %   GSCALE scales it to this range before processing.
     %
     %   The class of the output is the same as the class of the input.
     
     
     if length(varargin) == 0 % If only one argument it must be f.
        method = 'full8';
     else 
        method = varargin{1};
     end
     
     if strcmp(class(f), 'double') & (max(f(:)) > 1 | min(f(:)) < 0)
        f = mat2gray(f);
     end
     
     % Perform the specified scaling.
     switch method
     case 'full8'
        g = im2uint8(mat2gray(double(f)));
     case 'full16'
        g = im2uint16(mat2gray(double(f)));
     case 'minmax'
        low = varargin{2}; high = varargin{3};
        if low > 1 | low < 0 | high > 1 | high < 0
           error('Parameters low and high must be in the range [0, 1].')
        end
        if strcmp(class(f), 'double')
           low_in = min(f(:));
           high_in = max(f(:));
        elseif strcmp(class(f), 'uint8')
           low_in = double(min(f(:)))./255;
           high_in = double(max(f(:)))./255;
        elseif strcmp(class(f), 'uint16')
           low_in = double(min(f(:)))./65535;
           high_in = double(max(f(:)))./65535;    
        end
        % imadjust automatically matches the class of the input.
        g = imadjust(f, [low_in high_in], [low high]);   
     otherwise
        error('Unknown method.')
     end
    
  • 语法:
    g=gscale(f,method,low,high) g=gscale(f,method,low,high)

    • method取值:
      • ‘full8’(默认),省略low,high, 像素值标度到[0,255]
      • ‘full16’, 省略low,high, 像素值标度到[0,65535]
      • ‘minmax’, 必须有low,high ([0,1]之间)

直方图处理

  • 数字图像的直方图:
    h(rk)=nk h(r_k)=n_k

    • 离散函数
    • 灰度级范围 [0,G] ,共 L 个级别。
    • r_k 表示第 k 级灰度值,n_k 表示图像中灰度为 r_k 的像素个数。
  • 归一化直方图:
    p(rk)=nk/MN p(r_k)=n_k/MN

    • 它是各级灰度在图像中出现的频率统计
    • 若将灰度级 r_k 看成一个随机变量,则它是 r_k 的概率分布列的一个估计
  • 应用:

    • 图像增强
    • 图像压缩
    • 图像分割
  • 代码实现:
    M文件:

      function hist(f)
      %HIST gives the histogram of an image f
      
      % First convert the image to uint8
      f_uint8 = im2uint8(f);
      [M, N] = size(f);
      MN = M * N;
      r = 0:255;
      n = zeros(1,256);
      for i = 1:256
          n(i) = sum(sum(f_uint8 == (i-1)))/MN;
      end
      plot(r, n)
      xlim([0, 255])
    
  • IPT提供的API
    imhist:

     >> f = imread('Fig0306(a)(bone-scan-GE).tif');
     >> imhist(f)
    

    3.49

  • 例:
    3.21

    从上图中可以看出:

  • 暗图像的直方图集中在低灰度端

  • 亮图像的直方图集中在高灰度端

  • 低对比度图像的直方图集中在中间

  • 高对比度图像的直方图分布均匀,覆盖范围广

总结:
若一幅图像的像素分布接近于整个灰度级上的均匀分布,则图像会有高对比度,灰色调变化较大,最终效果是一幅灰度细节丰富,动态范围较大的图像。

直方图均衡

目标:找一个灰度变换,使得变换后的图像的直方图接近整个灰度范围上的均匀分布,以此得到高对比度,动态范围广的图像。

step1: 目标灰度变换需满足的必要条件:(假设灰度级归一化为 [0,1] )
s=T(r) s=T(r)

  • T(r)T(r) 在区间 [0,1] 上单调增
    这是为了防止出现灰度反转,造成人为缺陷。
  • T(r)T(r) 的值域也是 [0,1]
    这是保证输入输出的灰度范围相同

step2: 假设输入 r ,输出 s 均为 [0,1] 上随机变量,pr(r),ps(s)p_r(r),p_s(s) 分别是它们的概率密度函数,则令:
s=0rpr(w)dw s=\int_0^r p_r(w)dw
由概率密度函数非负和归一化条件知step1中条件满足

由概率论可知:
ps(s)=pr(r)drdsdsdr=pr(r)ps(s)=pr(r)1pr(r)=1 \begin{aligned} &amp; p_s(s)=p_r(r)|\frac{dr}{ds}|\\ &amp; \frac{ds}{dr}=p_r(r)\\ &amp; p_s(s)=p_r(r)\frac{1}{p_r(r)}=1 \end{aligned}
说明在上述灰度变换下,s服从 [0,1] 上的均匀分布

step3: 离散化

数字图像的灰度值总是离散的,用求和代替积分,求出离散形式的灰度变换函数
sk=T(rk)=j=1kpr(rj)=j=1knjn ,k=1,2,...L s_k=T(r_k)=\sum\limits_{j=1}^kp_r(r_j)=\sum\limits_{j=1}^k\frac{n_j}{n}\ ,k=1,2,...L
由于连续时, s 是 [0,1] 到 [0,1] 的单调增函数,离散化,再还原为原灰度范围(必要时四舍五入)后,sks_k 仍然满足条件

例:假设一幅大小为 64×6464\times 64 像素的3比特图像,灰度级为 [0,7] 中的整数,具体信息如下:

灰度级 归一化灰度级 rk nk pr(rk) sk sk(四舍五入)转化到原灰度级范围
0 0 r1 790 0.19 0.19 1.33取1
1 0.125 r2 1023 0.25 0.44 3.08取3
2 0.25 r3 850 0.21 0.65 4.55取5
3 0.375 r4 656 0.16 0.81 5.67取6
4 0.5 r5 329 0.08 0.89 6.23取6
5 0.625 r6 245 0.06 0.95 6.65取7
6 0.75 r7 122 0.03 0.98 6.86取7
7 0.875 r8 81 0.02 1 7.00取7

3.22

note: 离散化后,直方图只是概率密度的近似,不能证明离散直方图均衡能导出均匀的直方图,但这种变换有展开输入图像直方图的趋势,均衡化后的灰度级更宽了,最终效果是增强了对比度。

IPT提供的API:

histeq

 f = imread('Fig0316(4)(bottom_left).tif');
 g = histeq(f,256);
 figure
 subplot(2,2,1);imshow(f),title('原图');
 subplot(2,2,2);imhist(f),title('原图的直方图');
 subplot(2,2,3);imshow(g),title('均衡化后');
 subplot(2,2,4);imhist(g),title('均衡化后的直方图');

3.50

直方图匹配

目标:找一个灰度变换,使得变换后的图像的直方图是我们预定的直方图,这个预定是我们根据需求确定的。

step1: 先做直方图均衡,将r变换成s,s服从[0,1]上均匀分布
s=T(r) s=T(r)
step2: 假设z是服从我们预定的直方图的随机变量,则对z做直方图均衡,同样得到s
s=G(z) s=G(z)
step3: 求逆

要求G满足:

  • G(z) 在区间 [0,1] 上严格单调增
    为了保证G的反函数存在

  • G(z) 的值域也是 [0,1]
    这是保证输入输出的灰度范围相同

    求逆得:
    z=G1[T(r)] z = G^{-1}[T(r)]
    是我们要得变换

step4: 离散化

实践中我们处理的是离散值,最终得到预定直方图的一个近似

先有直方图均衡得到:
sk=T(rk)sk=G(zq) s_k=T(r_k)\\ s_k=G(z_q)
根据这两个对应关系可得:
zq=G1[T(rk)] z_q=G^{-1}[T(r_k)]
例:假设一幅大小为 64×6464\times 64 像素的3比特图像,灰度级为 [0,7] 中的整数,具体信息如下:

像素值 实际的pr 规定的pz r到s z到s r到z
0 0.19 0.00 1 0 3
1 0.25 0.00 3 0 4
2 0.21 0.00 5 0 5
3 0.16 0.15 6 1 6
4 0.08 0.20 6 2 6
5 0.06 0.30 7 5 7
6 0.03 0.20 7 6 7
7 0.02 0.15 7 7 7

实际中不需要G可逆,因为我们处理的是离散值,由r去找对应的s, 由s去找那些z与s对应,如果有多个,取与r最近的那个,如果没有相应的z, 就重新找与r对应的s, 找最接近的。
3.56

  • IPT提供的API:

    histeq:

     f = imread('Fig0323(a)(mars_moon_phobos).tif');
     ff = histeq(f,256);
     
     figure
     subplot(2,2,1);imshow(f),title('原图');
     subplot(2,2,2);imhist(f),title('原图的直方图');
     subplot(2,2,3);imshow(ff),title('均衡化后');
     subplot(2,2,4);imhist(ff),title('均衡化后的直方图');
    

    3.57

均衡化后效果并不好,因为原图集中在暗色区,使得离散化的灰度变换很快达到最大灰度值,均衡化后,集中到了亮色区

使用直方图匹配改进:

原图的直方图有两个峰值,大的靠近0,小的靠近255,我们用两个峰值的高斯函数拟合这种直方图图像

多峰值高斯函数:

 function p = twomodegauss(m1, sig1, m2, sig2, A1, A2, k)
 %TWOMODEGAUSS Generates a two-mode Gaussian function.
 %   P = TWOMODEGAUSS(M1, SIG1, M2, SIG2, A1, A2, K) generates a
 %   two-mode, Gaussian-like function in the interval [0,1].  P is a
 %   256-element vector normalized so that SUM(P) equals 1.  The mean
 %   and standard deviation of the modes are (M1, SIG1) and (M2,
 %   SIG2), respectively. A1 and A2 are the amplitude values of the
 %   two modes.  Since the output is normalized, only the relative
 %   values of A1 and A2 are important.  K is an offset value that
 %   raises the "floor" of the function.  A good set of values to try
 %   is M1=0.15, S1=0.05, M2=0.75, S2=0.05, A1=1, A2=0.07, and
 %   K=0.002.
 
 c1 = A1 * (1 / ((2 * pi) ^ 0.5) * sig1);
 k1 = 2 * (sig1 ^ 2);
 c2 = A2 * (1 / ((2 * pi) ^ 0.5) * sig2);
 k2 = 2 * (sig2 ^ 2);
 z  = linspace(0, 1, 256);

 p = k + c1 * exp(-((z - m1) .^ 2) ./ k1) + ...
     c2 * exp(-((z - m2) .^ 2) ./ k2);
 p = p ./ sum(p(:));
 p = twomodegauss(0.15, 0.05, 0.75, 0.05, 1, 0.07, 0.002);
 f = imread('Fig0323(a)(mars_moon_phobos).tif');
 g = histeq(f,p);
 figure
 subplot(2,2,1);imshow(f),title('原图');
 subplot(2,2,2);imhist(f),title('原图的直方图');
 subplot(2,2,3);imshow(g),title('均衡化后');
 subplot(2,2,4);imhist(g),title('均衡化后的直方图');

3.58

空间滤波基础

线性空间滤波器

在图像像素上执行线性操作的滤波器称为线性空间滤波器

一般形式:

假设领域大小取成 m×n=(2a+1)×(2b+1)m\times n=(2a+1)\times(2b+1) ,a,b 是正整数。
g(x,y)=u=au=av=bv=bw(u,v)f(x+u,y+v) g(x,y)=\sum\limits_{u=-a}^{u=a}\sum\limits_{v=-b}^{v=b}w(u,v)f(x+u,y+v)

示意图:

3.28

空间相关与卷积

术语:离散单位冲激:包含一个1,其余全是0的函数。

空间相关

相关操作的公式:

w(x,y)f(x,y)=u=au=av=bv=bw(u,v)f(x+u,y+v) w(x,y)☆f(x,y)=\sum\limits_{u=-a}^{u=a}\sum\limits_{v=-b}^{v=b}w(u,v)f(x+u,y+v)

  • 对所有位移量 x,y 求值
  • 在靠近边界处,需要填充 f
与离散单位冲击作用:
  • 离散单位冲击与滤波器如下:
    3.29
  • 进行相关操作:
    3.33
  • 裁剪后最终结果:
    3.34
    在1的位置生成一个翻转180度的滤波器

空间卷积

相关操作的公式:

w(x,y)f(x,y)=u=au=av=bv=bw(u,v)f(xu,yv) w(x,y)*f(x,y)=\sum\limits_{u=-a}^{u=a}\sum\limits_{v=-b}^{v=b}w(u,v)f(x-u,y-v)

  • 对所有位移量 x,y 求值
  • 在靠近边界处,需要填充 f
  • 相当于w先翻转180°,再做相关
与离散单位冲击作用:
  • 离散单位冲击与滤波器如下:
    3.29
  • 进行相关操作:
    3.33
  • 裁剪后最终结果:
    3.38
    在1的位置生成一个滤波器

代码实现

  function g = myImfilter(f, w, filter_mode, bd_type, size_option)
  %MYIMFILTER implements the linear space filtering operation of an image f
  % f is the input image
  % g is the output image
  % w is a linear space filter, which is given by users. The size of w must
  % be odd
  % filter_mode has two options, 'corr'and 'conv',which indicate 
  % the relevant and the convolution operation respectly
  % bd_type has four options
  % p , a value, means fill f with p(The default is 0)
  % 'rep' means fill f with the boundary value
  % 'sym' means image size expanding by mirroring their boundaries
  % 'cir' means fill f periodicly
  % size_option has two options
  % 'full' outputs extended g
  % 'same' outputs g with the same size of f
  
  % get the size oof w
  [m,n] = size(w);
  if mod(m*n,2)==0
      error('The size of w must be odd')
  end
  
  % choose flipping w or not based on filter_mode
  if strcmp(filter_mode,'corr')
      ww = w;
  else
      ww = rot90(w,2);
  end
  
  %fill f according to the option of bd_type
  a = (m-1)/2;
  b = (n-1)/2;
  [M,N] = size(f);
  ff = zeros(M + 2*a, N + 2*b);
  ff((a+1):(a+M),(b+1):(b+N)) = f(:,:);
  if strcmp(bd_type,'rep')
      for i = 1:a
          ff(i,(b+1):(b+N)) = f(1,:);
      end
      for i = (a+M+1):(2*a + M)
          ff(i,(b+1):(b+N)) = f(M,:);
      end
      for j = 1:b
          ff(:,j) = ff(:,b+1);
      end
      for j = (b+N+1):(2*b+N)
          ff(:,j) = ff(:,b+M);
      end
  elseif strcmp(bd_type,'sym')
      for i = 1:a
          ff(i,(b+1):(b+N)) = f(a+1-i,:);
      end
      for i = 1:a
          ff(a+M+i,(b+1):(b+N)) = f(M-i+1,:);
      end
      for j = 1:b
          ff(:,j) = ff(:,2*b+1-i);
      end
      for j = 1:b
          ff(:,b+N+j) = ff(:,b+N+1-j);
      end
  elseif strcmp(bd_type,'cir')
      for i = 1:a
          ff(i,(b+1):(b+N)) = f(M-a+i,:);
      end
      for i = 1:a
          ff(a+M+i,(b+1):(b+N)) = f(i,:);
      end
      for j = 1:b
          ff(:,j) = ff(:,N+j);
      end
      for j = 1:b
          ff(:,b+N+j) = ff(:,b+j);
      end
  else
      p = bd_type;
      for i = 1:a
          ff(i,(b+1):(b+N)) = ones(1,N).*p;
      end
      for i = 1:a
          ff(a+M+i,(b+1):(b+N)) = ones(1,N).*p;
      end
      for j = 1:b
          ff(:,j) = ones(M+2*a,1).*p;
      end
      for j = 1:b
          ff(:,b+N+j) = ones(M+2*a,1).*p;
      end
  end
  
  gg = zeros(M+2*a,N+2*b);
  for i = (a+1):(a+M)
      for j = (b+1):(b+N)
          gg(i,j) = sum(sum(ww.* ff((i-a):(i+a),(j-b):(j+b))));
      end
  end
      
  if strcmp(size_option, 'full')
      g = gg;
  else
      g = zeros(M,N);
      g = gg((a+1):(a+M),(b+1):(b+N));
  end

IPT提供的API

imfilter

  • 语法:
    g=imfilter(f,w,filter_mode,boundary_options,size_options) g = imfilter(f,w,filter\_mode,boundary\_options,size\_options)

    • f 是输入图像
    • w是滤波器
    • filter_mode 是选择操作模式,‘ corr ’ 代表相关,‘ conv ’ 代表卷积
    • boundary_options 是选择填充边界的方式,数值p表示用p填充,‘replicate’ 表示用边界值扩展,‘symmetric’ 表示镜像反射边界进行扩展,‘circular’ 表示将数组当成二维周期函数扩展
    • size_options 是选择输出模式,‘full’ 输出扩展的数组,‘same’ 输出与输入同样大小的裁剪数组
  • 不同模式处理结果对比:

     f = imread('Fig0216(a).tif');
     w = ones(31)./(31^2);
     g1 = imfilter(f,w,'corr','same');
     g2 = imfilter(f,w,'corr','replicate','same');
     g3 = imfilter(f,w,'corr','symmetric','same');
     g4 = imfilter(f,w,'corr','circular','same');
     figure
     subplot(3,2,[1,2]);imshow(f),title('原图');
     subplot(3,2,3);imshow(g1,[]),title('0填充');
     subplot(3,2,4);imshow(g2,[]),title('replicate');
     subplot(3,2,5);imshow(g3,[]),title('symmetric');
     subplot(3,2,6);imshow(g4,[]),title('circular');
    

    3.39
    此处w对称,所以相关和卷积结果一样,分析上图可以看出,w是一个平滑滤波器

    note:imfilter会将输出图像转化为输入图像的类型,比如,输入图像是‘uint8’,w没有归一化,计算可能超出255,最终会被imfliter截断,导致输出异常,如下图。因此,使用时注意将w归一化

     f = imread('Fig0216(a).tif');
     w = ones(31);
     g1 = imfilter(f,w,'corr','same');
     figure
     subplot(1,2,1);imshow(f),title('原图');
     subplot(1,2,2);imshow(g1,[]),title('0填充');
    

    3.40

平滑空间滤波器/均值滤波器/低通滤波器

  • 作用:模糊处理与降噪
  • w系数求和为1:滤波时,相当于对 (x,y) 邻域里的所有像素值做加权平均。
    特别的,w系数相同,相当于做简单平均,此时也称为盒状滤波器
  • 例:
     f = imread('Fig0333(a)(test_pattern_blurring_orig).tif');
     m = [3 5 9 15 35];
     figure
     subplot(2,3,1);imshow(f),title('原图');
     for i = 1:5
         w = ones(m(i))./(m(i)^2);
         g = imfilter(f,w,'corr','same');
         subplot(2,3,i+1);
         imshow(g),title(m(i));
     end
    
    3.51
    当图像有锯齿的时候,可以选择平滑滤波器,从上图可以看出,邻域取得越大,模糊效果越明显,需要根据图像情况,选择合适的滤波器大小。

锐化空间滤波器

作用:

突出图像中过渡部分

基础:

图像的一阶微分与二阶微分

一阶微分要求:

  • 在恒定灰度值区域,微分为零
  • 在灰度台阶或斜坡处非零
  • 沿着斜坡的微分值非零

二阶微分要求:

  • 在一阶微分值恒定区域,二阶微分为零
  • 在灰度台阶或斜坡起点处非零
  • 沿着斜坡二阶微分值为零

这些要求是从连续情形下,微分的性质类比而来,在离散情形下,此处的一阶微分实指差分

定义:

f 在 (x,y) 处的一阶微分:
fx=f(x+1,y)f(x,y)fy=f(x,y+1)f(x,y) \frac{\partial f}{\partial x}=f(x+1,y)-f(x,y)\quad \frac{\partial f}{\partial y}=f(x,y+1)-f(x,y)
f 在 (x,y) 处的二阶微分:
2f2x=f(x+1,y)+f(x1,y)2f(x,y)2f2y=f(x,y+1)+f(x,y1)2f(x,y) \frac{\partial^2 f}{\partial^2 x}=f(x+1,y)+f(x-1,y)-2f(x,y)\quad\frac{\partial^2 f}{\partial^2 y}=f(x,y+1)+f(x,y-1)-2f(x,y)

图像锐化的拉普拉斯算子

2f(x,y)=2f2x+2f2y=f(x+1,y)+f(x1,y)+f(x,y+1)+f(x,y1)4f(x,y) \nabla^2f(x,y) =\frac{\partial^2 f}{\partial^2 x}+\frac{\partial^2 f}{\partial^2 y}=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y)
术语:各向同性滤波器是指旋转不变的滤波器,w具有旋转不变性。

滤波器模板

由上述拉普拉斯算子可以得到如下滤波器模板
w=010141010 w= \begin{matrix} 0 &amp; 1 &amp; 0 \\ 1 &amp; -4 &amp; 1 \\ 0 &amp; 1 &amp; 0 \end{matrix}
这是一个以90°为增量的各向同性模板

例:找一个以45°为增量的各项同性模板

在拉普拉斯算子的基础上加入对角线方向,形成算子:
T=f(x+1,y+1)+f(x+1,y)+f(x1,y)+f(x1,y1)+f(x,y+1)+f(x,y1)+f(x+1,y1)+f(x1,y+1)8f(x,y) \begin{aligned} T &amp; =f(x+1,y+1)+f(x+1,y)+f(x-1,y)+f(x-1,y-1)\\ &amp; +f(x,y+1)+f(x,y-1)+f(x+1,y-1)+f(x-1,y+1)\\ &amp; -8f(x,y) \end{aligned}
模板为
w=111181111 w= \begin{matrix} 1 &amp; 1 &amp; 1 \\ 1 &amp; -8 &amp; 1 \\ 1 &amp; 1 &amp; 1 \end{matrix}

拉普拉斯图像锐化

原图像减去拉普拉斯变换后的图像可以复原背景特性并保持拉普拉斯锐化处理效果
g(x,y)=f(x,y)+c2f(x,y) g(x,y)=f(x,y)+c\nabla^2f(x,y)
若w中心系数为负,则取c=-1;反之,取c=1

例:

   f = imread('Fig0338(a)(blurry_moon).tif');
   w = [0 1 0;1 -4 1;0 1 0];
   delta_f = imfilter(f,w,'corr','same');
   delta_f_scale = gscale(delta_f,'full8');
   g = f-delta_f;
   figure
   subplot(1,4,1);imshow(f),title('月球北极的模糊图像');
   subplot(1,4,2);imshow(delta_f),title('未标定的拉普拉斯滤波后图像');
   subplot(1,4,3);imshow(delta_f_scale),title('标定的拉普拉斯滤波后图像');
   subplot(1,4,4);imshow(g),title('锐化后图像');

3.52

非锐化掩蔽与高提升滤波

  • 模糊原图像
  • 从原图像中减去模糊图像(称为模板)
  • 原图像加上权重模板
    公式:
    g(x,y)=f(x,y)+k(f(x,y)f(x,y)) g(x,y)=f(x,y)+k(f(x,y)-f^*(x,y))
  • f(x,y)f^*(x,y) 是模糊后的图像
  • k0k\ge 0
  • k=1时称为非锐化掩蔽;k>1时称为高提升滤波;k<1时不强调模板的贡献

图像说明:
3.53

术语:高斯平滑滤波器是指w由二维高斯函数生成的滤波器

例:

 f = imread('Fig0340(a)(dipxe_text).tif');
 % 生成高斯平滑滤波器,sigma=3
 w = zeros(5,5);
 for i = 1:5
     for j = 1:5
         w(i,j) = exp(-((i-3)^2+(j-3)^2)/18);
     end
 end
 ww = w./25;
 gauss_f = imfilter(f,ww,'corr','same')
 g_mask = f-gauss_f;
 g1 = f+g_mask;
 g45 = f+4.5.*g_mask;
 
 figure
 subplot(1,5,1);imshow(f),title('原图');
 subplot(1,5,2);imshow(gauss_f),title('高斯模糊后的结果');
 subplot(1,5,3);imshow(g_mask),title('模板');
 subplot(1,5,4);imshow(g1),title('非锐化掩蔽');
 subplot(1,5,5);imshow(g45),title('k=4.5的高提升滤波');

3.54

梯度法

  • 图像梯度:
    f=[fx,fy]T \nabla f=[\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}]^T

  • 梯度幅值:
    M(x,y)=(fx2+fy2) M(x,y)=\sqrt{(f_x^2+f_y^2)}
    图像 M(x,y) 称为梯度图像

  • 例:
    sobel算子:(3×33\times 3模板)

    • 梯度近似:
      gx=(z7+2z8+z9)(z1+2z2+z3)gy=(z3+2z6+z9)(z1+2z4+z7) g_x=(z_7+2z_8+z_9)-(z1+2z_2+z_3)\\ g_y=(z_3+2z_6+z_9)-(z1+2z_4+z_7)
      • x方向用第3行减第1行近似梯度
      • y方向用第3列减第1列近似梯度
      • 中间系数用2是为了保证模板系数和为0,从而满足灰度恒定区域,梯度为0
    • 模板
      wx=[121000121] w_x= \begin{bmatrix} -1 &amp; -2 &amp; -1 \\ 0 &amp; 0 &amp; 0 \\ 1 &amp; 2 &amp; 1 \end{bmatrix}
      wy=[101202101] w_y= \begin{bmatrix} -1 &amp; 0 &amp; 1 \\ -2 &amp; 0 &amp; 2 \\ -1 &amp; 0 &amp; 1 \end{bmatrix}
      这两个模板称为soble算子
  • 梯度图像
    M(x,y)(z7+2z8+z9)(z1+2z2+z3)+(z3+2z6+z9)(z1+2z4+z7) M(x,y) \approx |(z_7+2z_8+z_9)-(z1+2z_2+z_3)|+|(z_3+2z_6+z_9)-(z1+2z_4+z_7)|

  • 例:使用梯度进行边缘增强

       f = imread('Fig0342(a)(contact_lens_original).tif');
       % soble算子
       wx = [-1 -2 -1;0 0 0 ;1 2 1]./9;
       wy = [-1 0 1;-2 0 2;-1 0 1]./9;
       fx = imfilter(f,wx,'corr','same');
       fy = imfilter(f,wy,'corr','same');
       M = abs(fx)+abs(fy);
       
       figure
       subplot(1,2,1);imshow(f),title('原图');
       subplot(1,2,2);imshow(M),title('梯度图像');
    

    3.55
    原图是一个玻璃镜片,右下角有一个缺口,进行边缘增强后,缺口更明显

非线性空间滤波器

统计排序滤波器

  • 中值滤波器:(x,y) 处的输出像素值是它邻域内像素值的中值
  • 最大值滤波器: (x,y) 处的输出像素值是它邻域内像素值的最大值
  • 最小值滤波器:(x,y) 处的输出像素值是它邻域内像素值的最小值

中值滤波器有很好的去噪能力

术语:椒盐噪声:以黑白点的形式叠加在图像上

IPT的标准空间滤波器

fspecial

生成线性滤波器

  • 语法:
    w=fspecial(type,parameters) w = fspecial(&#x27;type&#x27;,parameters)
    参数取值如下:
    3.59
    3.60

  • 用fspecial生成一个拉普拉斯滤波器,对图像进行锐化
     w = fspecial('laplacian',0);
     f = imread('Fig0338(a)(blurry_moon).tif');
     ff = im2double(f);
     g = imfilter(ff,w,'replicate');
     g2 = ff-g;
     figure
     subplot(1,3,1);imshow(f),title('月球北极的模糊图像');
     subplot(1,3,2);imshow(g),title('拉普拉斯滤波后图像');
     subplot(1,3,3);imshow(g2),title('图像锐化'); 
    
    3.61

ordfilt2

生成统计排序滤波器

  • 语法:
    g=ordfilt2(f,ord,domain) g = ordfilt2(f,ord,domain)
    • f输入图像
    • ord取排序后第几百分位代替邻域中心处的值
    • domain 由0,1组成的 m×nm\times n 矩阵,0处的值不参与排序
  • 例:
    • 最大值滤波器:
      g=ordfilt2(f,mn,ones(m,n)) g = ordfilt2(f,m*n,ones(m,n))
    • 最小值滤波器
      g=ordfilt2(f,1,ones(m,n)) g = ordfilt2(f,1,ones(m,n))
    • 中值滤波器
      g=ordfilt2(f,median(1:mn),ones(m,n)) g = ordfilt2(f,median(1:m*n),ones(m,n))

medfilt2

二维中值滤波函数

  • 语法:
    g=mefilt2(f,[m,n],padopt) g = mefilt2(f,[m,n],padopt)
    • [m,n]邻域大小
    • padopt
      • ‘zeros’ 零填充(默认)
      • ‘symmetric’ 镜像对称填充
      • ‘indexed’ 如果是double类,则1填充,否则0填充
  • 例:去除椒盐噪声
     f = imread('Fig0335(a)(ckt_board_saltpep_prob_pt05).tif');
     g = medfilt2(f);
     g2 = medfilt2(f,'symmetric');
     
     figure
     subplot(1,3,1);imshow(f),title('原图');
     subplot(1,3,2);imshow(g),title('0填充');
     subplot(1,3,3);imshow(g2),title('镜像对称填充');
    
    3.62

参考文献:

《数字图像处理》第三版

《数字图像处理MATLAB版》

2012-01-09 09:40:40 iteye_10132 阅读数 40
  • OpenCV3.2 Java图像处理视频学习教程

    OpenCV3.2 Java图像处理视频培训课程:基于OpenCV新版本3.2.0详细讲述Java OpenCV图像处理部分内容,包括Mat对象使用、图像读写、 基于常用核心API讲述基本原理、使用方法、参数、代码演示、图像处理思路与流程讲授。主要内容包括opencv像素操作、滤波、边缘提取、直线与圆检测、形态学操作与分水岭、图像金子塔融合重建、多尺度模板匹配、opencv人脸检测、OpenCV跟Tomcat使用实现服务器端图像处理服务。

    4137 人正在学习 去看看 贾志刚
图像滤波
2010年09月29日
  刚获得的图像有很多噪音。这主要由于平时的工作和环境引起的,图像增强是减弱噪音,增强对比度。
  想得到比较干净清晰的图像并不是容易的事情。为这个目标而为处理图像所涉及的操作是设计一个适合、匹配的滤波器和恰当的阈值。
  常用的有
  高斯滤波
  均值滤波
  中值滤波
  最小均方差滤波
  Gabor滤波
  由于[b]高斯函数的[/b]傅立叶变换仍是高斯函数, 因此高斯函数能构成一个在频域具有平滑性能的低通滤波器。可以通过在频域做乘积来实现高斯滤波。
  [b]均值滤波[/b]是对是对信号进行局部平均, 以平均值来代表该像素点的灰度值。
  矩形滤波器(Averaging Box Filter)对这个二维矢量的每一个分量进行独立的平滑处理。通过计算和转化 ,得到一幅单位矢量图。
  这个 512×512的矢量图被划分成一个 8×8的小区域 ,再在每一个小区域中 ,统计这个区域内的主要方向 ,亦即将对该区域内点方向数进行统计,最多的方向作为区域的主方向。
  于是就得到了一个新的64×64的矢量图。这个新的矢量图还可以采用一个 3×3模板进行进一步的平滑。
  中值滤波是常用的非线性滤波方法 ,也是图像处理技术中最常用的预处理技术。
  它在平滑脉冲噪声方面非常有效,同时它可以保护图像尖锐的边缘。加权中值滤波能够改进中值滤波的边缘信号保持效果。但对方向性很强的指纹图像进行滤波处理时 ,有必要引入方向信息,即利用指纹方向图来指导中值滤波的进行。
  最小均方差滤波器,亦称维纳滤波器,其设计思想是使输入信号乘响应后的输出,与期望输出的均方误差为最小。
  Gabor变换是英国物理学家 Gabor提出来的,由“测不准原理”可知,它具有最小的时频窗,即Gabor函数能做到具有最精确的时间-频率的局部化;另外, Gabor函数与哺乳动物的视觉感受野相当吻合,这一点对研究图像特征检测或空间频率滤波非常有用。恰当的选择其参数, Gabor变换可以出色地进行图像分割、识别与理解。如文献提出的基于Gabor滤波器的增强算法。
  评论这张
  
  
  
  转发至微博

滤波的分享

阅读数 120

没有更多推荐了,返回首页