看到太多的人把MATLAB当成C语言来写,甚至像矩阵求和这种的,都要写个for循环来做,不仅整个程序看起来很乱,而且不管开发效率还是程序执行效率都降低了太多。我在研一学人工神经网络这门课程的时候,就有看过别人写的一些MATLAB程序,用矩阵的一些操作代替了for循环,程序看起来简洁干净,效率也有提高。
这篇文章接下来会介绍几个函数,能把原来需要用for循环来处理矩阵的,通过自带函数替换掉for循环。用这些函数的好处是代码会简洁很多,同时MATLAB也会优化操作过程,提高效率。

1. 对矩阵非线性滤波操作或滑动窗口操作

1.1 colfilt

函数语法:
     B = colfilt(A,[m n],block_type,fun)
     B = colfilt(A,[m n],[mblock nblock],block_type,fun)
     B = colfilt(A,’indexed’,…)

其中:

  • A - 要进行滤波的矩阵
  • [m, n] - 窗口大小
  • block_type:
    • ‘sliding’ - 滑动窗口模式,逐个像素滑动,下面主要讲的是这个模式
    • ‘distinct’ - 块模式,对矩阵使用fun函数按指定窗口大小操作
  • fun - 滤波函数
  • parameters - 函数要传入的其他参数

函数功能:
     对矩阵的领域操作

当block_type为‘sliding’模式时,colfilt函数执行过程如下:

  1. 首先对f进行列变换,假设f是MxN的图像,则会生成一个 mn x MN 的矩阵。其实就是把原图像每个像素邻域内的mn个像素排成一列,原图像总共有MN个像素,因此变换后的矩阵就有MN

  2. 然后colfilt会每次选择多列作为一个矩阵传入到fun中进行处理,到底选择几列是由colfilt按照一定策略进行的。而传入矩阵的行数肯定是mn。

  3. fun函数对传入的矩阵进行运算,一般是一列一列的处理,然后返回一个行向量,这个行向量的每个元素对应每列处理的结果,因此传入的矩阵有多少列就按序生成多少个结果放到对>

  4. colfilt对返回的行向量进行存储、反变换,最后得到和原图像等大小的一个矩阵,这就是非线性滤波的最终结果。

示例,对图像进行中值滤波

1
2
3
4
5
I = imread('tire.tif');
I2 = uint8(colfilt(I,[5 5],'sliding',@median));
figure
subplot(1,2,1), imshow(I), title('Original Image')
subplot(1,2,2), imshow(I2), title('Filtered Image')

1.2 nlfilter

函数语法:
     B = nlfilter(A, [m n], fun)
     B = nlfilter(A, ‘indexed’,…)
其中:

  • A - 要进行滤波的矩阵
  • [m, n] - 窗口大小
  • fun - 滤波函数

函数功能:
     使用滑动窗口操作矩阵

nlfilter比colfilt占用更少的内存,但执行效率没有nfilter好。图像处理过程中,使用滑动窗口操作的过程都可可以使用这两个函数。如,对图像应用滑动窗口技术,求解每个窗口内>

1
2
Iout = nlfilter(image, [8 8], @std2);
Iout = colfilt(image, [8 8], 'sliding', @std);

2. 块处理函数 blockproc

函数语法:
     B = blockproc(A,blockSize,fun)
     B = blockproc(src_filename,blockSize,fun)
     B = blockproc(adapter,blockSize,fun)
     blockproc(___,Name,Value,…)
其中:

  • A - 要进行块操作的矩阵
  • blockSize - 块大小
  • fun - 处理函数

函数功能:
     将矩阵A分成blockSize大小的离散小块,然后对每一个小块使用fun函数进行处理

应用示例:

示例1:计算图像的缩略图

1
2
3
4
5
6
7
fun = @(block_struct) imresize(block_struct.data,0.15);
I = imread('pears.png');
I2 = blockproc(I,[100 100],fun);
figure;
imshow(I);
figure;
imshow(I2);

示例2:将RGB图像的红色和绿色通道互换

1
2
3
4
5
6
7
8
9
10
I = imread('peppers.png');
fun = @(block_struct) block_struct.data(:,:,[2 1 3]);
blockproc(I,[200 200],fun,'Destination','grb_peppers.tif');
figure;
imshow('peppers.png');
figure;
imshow('grb_peppers.tif');

% 这个其实可以直接转换
grb_peppers = I(:,:,[2 1 3]);

示例3:将图像缩小为原来的1/2,缩小后的图像像素为对应原始图像中4个像素值的平均

1
2
3
4
5
6
7
fun = @(x) mean(x.data(:));
I = imread('moon.tif');
I2 = blockproc(I,[2 2],fun);
figure;
imshow(I);
figure;
imshow(I2, []);

3. 向量化编程的利器 bsxfun

函数语法:
     C = bsxfun(fun, A, B)

函数功能:
     将两个矩阵中对应位置的元素应用fun函数处理,fun可以是系统内建的函数,注意:不能使用符号如”+”、”-“等,而应使用对应的函数句柄”@plus”、”@minus”,具体参看MATLAB帮助doc bsxfun

如果A、B两个矩阵一样大,可以直接使用A+B的方式来实现,但当A或B中有任何一维的长度是1时,函数会自动将该为扩展为和另一个矩阵对应维度的长度,并且,这样的扩展是虚拟的,会比使用repmat效率要高,占用内存也少。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
>> a = [1 2 3]

a =

1 2 3

>> b = [2; 4; 6]

b =

2
4
6

>> c = bsxfun(@plus, a, b)

c =

3 4 5
5 6 7
7 8 9

在解决下面的示例问题中,对比使用for循环与bsxfun运行效率:
有A,B两个矩阵,在这两个矩阵中,每行都代表一个样本,每列都表示一个特征,计算其高斯核
$$k(||x-xc||) = exp(-\frac{||x - xc||^2}{2\sigma^2})$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
A = rand(2000, 1000);
B = rand(3000, 1000);
beta = 0.5;
tic
% 使用for版本
K1 = zeros(size(A,1), size(B,1));
for i = 1 : size(A, 1)
for j = 1 : size(B, 1)
K1(i, j) = exp(-sum((A(i,:)-B(j,:)).^2)/beta);
end
end
toc

% 使用bsxfun版本
tic;
sA = (sum(A.^2, 2));
sB = (sum(B.^2, 2));
K3 = exp(bsxfun(@minus,bsxfun(@minus,2*A*B', sA), sB')/beta);
toc;

在我的电脑上,运行上述程序for版本耗时185.981322s,bsxfun版本的耗时0.542689s,bsxfun版本比for版本效率提升了超过300倍。

插曲:在实现bsxfun版本时,我最开始是用下面的这种方法

1
2
rowNum = size(A, 2);
K2 = exp(-sum(bsxfun(@minus, reshape(A, [], 1, rowNum), reshape(B, 1, [], rowNum)).^2, 3) / beta);

最后还出现内存不足,电脑卡死的状况。后来分析发现,我的这种方法是先把2维矩阵扩展成3维对应相减,然后再把结果求平方再求和,虽说bsxfun会在扩展矩阵的时候做优化,但结束之后返回的三维矩阵也非常非常大,耗尽内存。

4. arrayfun

函数语法:
     [B1,…,Bm] = arrayfun(func,A1,…,An)
     [B1,…,Bm] = arrayfun(func,A1,…,An,Name,Value)

函数功能:
     使用函数处理数组中每个元素,可以传入多个矩阵,则将多个矩阵同一位置的元素作为参数传入函数fun中处理。

示例程序:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
>> s(1).f1 = rand(3, 6);
>> s(2).f1 = magic(12);
>> s(3).f1 = ones(5, 10);
>> counts = arrayfun(@(x) numel(x.f1), s) % 计算每个元素的f1属性中元素个数

counts =

18 144 50

>> [nrows, ncols] = arrayfun(@(x) size(x.f1), s)

nrows =

3 12 5


ncols =

6 12 10
A =

0.1662 0.7518 0.2118
0.0817 0.5420 0.9325
0.9440 0.7209 0.7581


B =

0.0392 0.2084 0.6986
0.1176 0.4304 0.7956
0.4965 0.9252 0.3083

>> counts = arrayfun(@(x, y) x + y, A, B)

counts =

0.2054 0.9602 0.9104
0.1993 0.9724 1.7280
1.4405 1.6461 1.0664

5. cellfun

函数语法:
     [A1,…,Am] = cellfun(func,C1,…,Cn)
     [A1,…,Am] = cellfun(func,C1,…,Cn,Name,Value)

函数功能:
     使用函数处理cell中的数据。

cellfun函数的用法与arrayfun函数的用法类似,除了前者是用于处理cell,后者是用于处理矩阵。
示例程序参考MATLAB的帮助doc cellfun

6. spfun

函数语法:
     f = spfun(fun,S)

函数功能:
     对稀疏矩阵中的所有非零元素调用函数处理。

示例程序参考MATLAB的帮助doc spfun

7. structfun

函数语法:
     [A1,…,An] = structfun(func,S)
     [A1,…,An] = structfun(func,S,Name,Value)

函数功能:
     对结构体重所有元素调用函数处理。

示例程序参考MATLAB的帮助doc structfun