第一章 绪言

数字图像处理

数字图像处理(Digital Image Processing)是通过计算机对图像进行去除噪声、增强、复原、分割、提取特征等处理的方法和技术。

MATLAB

MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。

数字图像表示

一幅图像可以定义成一个二维函数f(x,y)，其中x和y是空间(平面)坐标，且任何坐标对(x,y)处的幅度f成为图像在这一点的亮度。术语灰度级通常指单色图像的亮度。彩色图像是由单色图像组合而成的。

坐标约定

采样和量化得到的是一个实数矩阵。我们对一幅图像f(x,y)采样后，可以得到一个M行、N列的图像。其中坐标的值是离散量。为了表示方便，坐标取整数值。

图像处理工具箱中用(r,c)表示行与列。同时，坐标的原点是(1,1)。因此，r的范围是1到M，c的范围是1到N，以整数递增。

作为矩阵的图像

通常把图像用M×N矩阵表示，M为行，N为列。矩阵的每个元素称为像素。

图像的输入/输出和显示

• imread(‘filename‘)
  输入函数，其中filename是一个包含文件全名的字符串。例如

  >> f = imread(‘chestxray.jpg‘);

• imshow(f)
  输出函数，其中f是图像数组。例如

>> f = imread(‘rose_512.tif‘);
>> imshow(f)

• 如果想在保留第一幅图像并输出第二幅图像，可以使用函数 figure ：

>> figure,imshow(g)

• imwrite(f,‘filename‘,‘quality‘,q)
  输出函数，这个函数可以将图像写入当前目录。例如

>> imwrite(g,‘filename.tif‘,‘quality‘,0);

类和图像类型

其中uint8和logical广泛地用于图像处理。

工具箱支持4种图像类型：

• 灰度级图像
• 二值图像
• 索引图象
• RGB图像

大多数单色图像的处理运算都是使用二值图像或灰度级图像来进行的。

灰度级图像

一幅灰度级图像是一个数据矩阵，矩阵的值表示灰度的浓淡。

二值图像

一幅二值图像是一个取值只有0和1的逻辑数组。其中0是黑色，1是白色。其他类型的图像也可以通过类型转换转换成二值图像。
如果A是一个由0和1构成的数组，那么可以如下创建一个逻辑数组B：
B = logical(A)
还有一些其他相关函数

• islogical(C)
  检测一个数组是否是logical类
• B = class_name(A)
  通用的类转换语法。class_name可以是im2uint8、im2uint16、im2double、im2single或mat2gray。

M函数编程

M文件

MATLAB中的M文件可以是简单地执行一系列MATLAB语句的脚本，也可以是能够接受参量产生一个或多个输出的函数。

M文件的定义如下：

• 函数定义行
• H1行
• 帮助文本
• 函数体
• 注释

函数定义行具有如下形式：

function [outputs] = name(inputs)

H1行是单独注释行，记录了概要信息。
帮助文本是紧跟在H1行后面的文本块，两者之间没有空行。帮助文本用来为函数提供注释和在线帮助。
函数体包含了执行计算和给输出变量赋值的所有MATLAB代码。
符号“%”后面的非H1行或帮助文本的所有行，都被认为是函数注释行。

算数运算符

使用方法与其他编程语言类似。

关系运算符

逻辑运算符

流控制

包括分支与循环等8个流控制语句。

数组索引

与其他语言不同的是，MATLAB中数组的下标是从1开始计算的。
一个1×N的数组称为行向量，一个M×1的数组称为列向量。例如：

>> v = [1 3 5 7 9]
v =
    1 3 5 7 9
>> v(2)
ans = 
    3

使用转置运算符(.‘),可将行向量转换为列向量(反之亦然)：

>> w = v.‘
w =
    1
    3
    5
    7
    9

要访问元素快，可使用MATLAB的冒号。例如：

>> v(1:3)
ans =
    1 3 5
>> v(3:end)
ans = 
    5 7 9
>> v([1 4 5])
ans = 
    1 7 9
>> v(1:2:end)
ans =
    1 5 9

类似地，矩阵的定义方法为

>> A = [1 2 3; 4 5 6; 7 8 9]
A = 
    1 2 3
    4 5 6
    7 8 9

同时，我们可以按处理向量的方式来选取矩阵中的元素，但是需要两个索引：一个对应于行，一个对应于列。也可以私用一个冒号作为索引来选择整行、整列或整个矩阵。

>> A(2,:)

ans =

     4     5     6

>> sum(A(:))

ans =

    45

>> A(:,2)

ans =

     2
     5
     8

另一种相当有用的索引形式是逻辑索引(logical indexing)。逻辑索引具有形式A(D)，其中A是一个数组，D是与A相同大小的逻辑数组。表达式A(D)提取A中与D中的1值元素相对应的所有元素。例如：

>> D = logical([1 0 0;0 0 1; 0 0 0])

D =

  3×3 logical 数组

   1   0   0
   0   0   1
   0   0   0

>> A(D)

ans =

     1
     6

对图像处理有用的最后一类索引是线性索引(linear indexing)。线性索引表达式使用单一下标来索引矩阵或高维数组。对于一个M×N的矩阵，元素(r,c)可用单一下标r+M(c-1)来访问。

函数句柄

命名函数句柄

例如：

>> f = @sin

f =

  包含以下值的 function_handle:

    @sin

>> f(pi/4)

ans =

    0.7071

>> sin(pi/4)

ans =

    0.7071

匿名函数句柄

通用格式为

@(input-argument-list) expression

例如：

>> g = @(x) x.^2

g =

  包含以下值的 function_handle:

    @(x)x.^2

>> g(2)

ans =

     4

代码优化

原始函数为 sinfun1.m
这个函数用来计算 f(x)=sin(x/(100*pi))

function [ y ] = sinfun1( M )
x = 0:M -1;
for k = 1 : numel(x)
    y(k) = sin(x(k)/(100*pi));
end

进行测速

>> M = 100;
>> f = @() sinfun1(M);
>> timeit(f)

ans =

   4.7852e-05

这是函数在M=100时执行所需的时间。继续测量M=500，1000，1500，……，20000时的执行时间：

>> M = 500:500:20000;
>> for k = 1:numel(M)
f = @() sinfun1(M(k));
t(k) = timeit(f);
end

>> t

t =

  1 至 9 列

    0.0001    0.0002    0.0002    0.0004    0.0004    0.0005    0.0004    0.0005    0.0007

  10 至 18 列

    0.0008    0.0006    0.0006    0.0007    0.0007    0.0008    0.0009    0.0009    0.0013

  19 至 27 列

    0.0010    0.0010    0.0011    0.0014    0.0012    0.0012    0.0013    0.0013    0.0014

  28 至 36 列

    0.0019    0.0020    0.0015    0.0020    0.0019    0.0017    0.0017    0.0020    0.0018

  37 至 40 列

    0.0019    0.0020    0.0020    0.0020

在这里，这段代码中，我们预测时间应该是线性的增长。因为与计算的M值有关。但是，在MATLAB中实际运行时，这个时间几乎是成M^2增长的。这是因为在函数中，每进行一次循环，y数组的大小都进行了隐式数组增长，MATLAB必须重新为其分配内存空间，这个过程中，重新分配和复制的开销很大，远超过计算sin值本身的时间。所以，我们在这个函数的基础上进行了一定的优化。
优化后新的函数如下：

function [ y ] = sinfun2( M )
x = 0:M -1;
y = zeros(1,numel(x));
for k = 1 : numel(x)
    y(k) = sin(x(k)/(100*pi));
end

这个函数中，使用了预分配。比较两次的所用的时间：

>> timeit(@() sinfun1(20000))

ans =

    0.0029

>> timeit(@() sinfun2(20000))

ans =

   5.3685e-04

可以看到运行时间的差别。
同样的，我们也可以用向量化方法来编写这个函数。

function [ y ] = sinfun3( M )
x = 0:M -1;
y = sin(x ./ (100*pi));
end

我们对其进行测速比较：

>> timeit(@() sinfun2(20000))

ans =

   6.5941e-04

>> timeit(@() sinfun3(20000))

ans =

   4.3185e-04

>> 

可以看出，sinfun2和sinfun3函数的速度相差不大，但都远远快于sinfun1函数。所以，不带循环的向量化方法同样可以加快运行速度。

向量化说明和函数meshgrid的介绍

下面，用一个MATLAB函数的两个版本，基于如下公式创建一幅图像：
f(x,y) = Asin(u0x+v0y)
第一个函数，用两层嵌套的for循环实现：

function f = twodsin1( A, u0, v0, M, N )
f = zeros(M,N);
for c = 1:N
    v0y = v0 * (c - 1);
    for r = 1:M
        u0x = u0 * (r - 1);
        f(r,c) = A*sin(u0x + v0y);
    end
end
end

同时，测量其运行时间

>> timeit(@() twodsin1(1,1/(4*pi),1/(4*pi),512,512))

ans =

    0.0091
    
>> timeit(@() twodsin0(1,1/(4*pi),1/(4*pi),512,512))

ans =

    0.0131

与不预先分配空间相比，显然提升了速度。
使用以下语句生成并显示图像

>> f = twodsin1(1,1/(4*pi),1/(4*pi),512,512);
>> imshow(f,[])

可以得到如下正弦图像

然后，使用meshgrid函数重写这个函数。这个函数的语法为

[C, R] = meshgrid(c, r)

函数meshgrid把坐标向量变换为两个数组C和R，这两个数组可以用来计算有两个变量的函数。
改写代码如下：

function f = twodsin2( A, u0, v0, M, N )
r = 0:M - 1;
c = 0:N - 1;
[C, R] = meshgrid(c, r);
f = A*sin(u0*R + v0*C);
end

测量运行时间

>> timeit(@() twodsin2(1,1/(4*pi),1/(4*pi),512,512))

ans =

    0.0051

相比之前，可以看到，向量化形式运行时间减少了50%以上。