Matlab矢量化编程技巧
?
级别：★经常看到一些Matlab初学者写出带有层层循环的代码，这些代码往往运行得很慢，而且非常难懂。Matlab提供了大量的命令来避免循环建议是：在确定要写一个针对矩阵操作的2重循环之前，请仔细阅读帮助中的Maximizing MATLAB Performance一节。里面有很多例子，这里举一个我自己碰到的问题：要统计数字图像的亮度最大值和最小值，由于数字图像可能是一个二维或三维的矩阵，开始我写出了如下的代码：switch ndims(img)case 2m = max(max(img));case 3m = max(max(max(img)));end后来对这种做法很不满意，如果有一个8维的矩阵，难道我要写max(max(max(max(max(max(max(max(x))))))))才行？后来发现无论矩阵是什么维数，实际上max( x(:) )就够了。点评：y = x(:) 这种方式可以返回矩阵的所有元素，形成一个列向量 y，无视 x 的维数。注：y 对 x 以列优先的方式排列 级别：★★★问题：提取一幅RGB图象中的蓝色部分解答：[Copy to clipboard] [ - ]CODE:A=imread('input.bmp');H = size(A);BB=repmat(255, H(1)*H(2), H(3));C=(A(:,:,1)==0 & A(:,:,2)==0 & A(:,:,3)==255);BB(C,1) = 0;BB(C,2) = 0;B = reshape(BB, H(1), H(2), H(3));imwrite(B,'output.bmp');点评：注意 repmat 和 reshape 的用法。解答的关键是对 reshape 命令的灵活使用，把三维矩阵转换为二维矩阵，从而使得求出满足条件的下标矩阵 C 可以对原矩阵 A 进行引用。否则，这类问题一定摆脱不了循环语句。 级别：★问题：赋值的技巧——当矩阵很多元素是相同时（下面以0为例，且不考虑稀疏矩阵的存储方式），“先生成全零阵，后再对特殊元素赋值”的方式将有利于缩短运行的时间。例如：把矩阵 A 作 0 扩展（变身：1个变4个，原来真身维持不变，其他假身为0），l例如：A = 1 2 3 4 -1 -2 -3 -4 0 4 2 9希望得到：B= 1 2 3 4 0 0 0 0 -1 -2 -3 -4 0 0 0 0 0 4 2 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0解答：[Copy to clipboard] [ - ]CODE:B = repmat(zeros(size(A)), 2, 2);B(1:size(A,1), 1:size(A,2)) = A;不推荐的做法：[Copy to clipboard] [ - ]CODE:[row, col] = size(A);for i = 1:2*row for j = 1:2*col if i >=1 && i <= row && j >= 1 && j <= col B(i, j) = A(i, j); else B(i, j) = 0; endend点评：也许是受制于 C 语言的影响，初学者往往把 matlab 当作 C 语言来使用，加上 matlab 宣称变量不需要定义即可以直接使用，因此导致大量的循环+单变量判断语句。掌握“：”运算符是向 matlab 矢量化编码进军的必经之路。 级别：★★问题：有个80*80的矩阵，如何去掉全为0的行或者列解答：[Copy to clipboard] [ - ]CODE:a = 0 0 0 -3 0 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0>> a(:,~sum(abs(a),1))=[]a = 0 0 -3 1 0 0 0 2 0 0 0 0 0 0 0>> a(~sum(abs(a),2),:)=[]a = 0 0 -3 1 0 0 0 2 0点评：又是一个摆脱循环思想的典型例子。除了“: ”运算符的妙用之外，逻辑法对矩阵的下标索引也相当重要。此外，线性代数的知识也在这里得到体现：|a|+|b|+|c|=0 当且仅当 a = b = c = 0 级别：★★★★★问题：对一幅灰度图像求它的二维联合分布密度——假设不使用 matlab 自带的函数直接求出来（如果有的话）问题详细描述：QUOTE:设 (i,j), (k,l) 两个任意象素点上的灰度值分别为 g(i,j) 和 g(k,l)，其值分别假设为 a, b，则图像灰度值的联合分布密度表示为P(a,b)=Pk{g(i,j)=a, g(k,l)=b}，其中 a, b 均为 0 到 L-1 之间的灰度整量电平。由于灰度图象是 unint8 型，那它的灰度级就是 256 个灰度级。根据上面关于二维联合分布密度的定义，即是要求得到这幅图中同时满足某个给定条件的点的数量，这个条件是两个位置上像素的灰度值与当前统计的灰度配对是一致的。整个二维联合分布密度矩阵与协方差矩阵类似。常规解答：[m,n]=size(a);a1=0for ii=1:256 for jj=1:256 I0(ii,jj)=0; for i=1:m for j=1:n for k=1:m for l=1:n if (a(i,j)==(ii-1))&(a(k,l)==(jj-1)) I0(ii,jj)=I0(ii,jj)+1; a1=a1+1; end end end end end endend;简短说明：程序中的 ii, jj 即是相关的 a, b 值，它可以为 0 到 255 之间的任意数（也可以理解为用 ii 作为横坐标， jj 作为纵坐标）。设置了 a, b 的值后，就用个 FOR 循环来寻找图象中任两个位置上能够同时满足像素值为 a 和像素值为 b 的点（即为 IF 语句后的条件），当满足条件时，那么在该点上加一，就能统计出同时满足 a, b 的点的数量。因为 a, b 可以取不同的值，那么就能得出一个 256*256 大小的矩阵，也就是想要的二维联合分布密度。优化解答：A = [1,2,3;2,3,6;6,4,5;1,1,3];B = unique(A);C = hist(double(A(:)), double(length(B)));D = C'*C;E = -diag(diag(D))+D;E = diag(C) + E;AE点评：这是 matlab 编程思想转变的一个典型例子，即从算法上去改进程序，使之减少循环。常规解法使用了6个循环，而优化解法则没有使用循环。转化的思想是：考虑到问题本身的特点，即：灰度配对的点的个数（简称点数） = 各灰度点数之积，因此首先求出矩阵 A 中互异的元素 B，然后，利用 hist 函数统计每个互异元素出现的次数 C，即各灰度点数。在此基础上，C 的转置与 C 的乘积矩阵 D 就是两互异灰度的点数。剩下的问题是把 D 的主对角线换掉（因为相乘时统计多了），这通过先置零，后添加来进行。 级别：★问题：求出满足特定条件的元素，并按行统计个数解答：不推荐的做法：[Copy to clipboard] [ - ]CODE:[m,n]=size(t2);for j=1:n for i=1:m if (t2(i,j)-y2(i,j))/t2(i,j)<0.2 hege(j)=hege(j)+1; end endend推荐的做法：[Copy to clipboard] [ - ]CODE:A = ((t2-y2)./t2) < 0.2;hege = sum(A,1);点评：充分发挥 matlab 基于矩阵的思想，从整体出发考虑问题，避免了使用基于 C 语言 的循环语句。 级别：★★★★★问题：游程平滑算法：将连续的且个数小于某个阈值的0全部替换成1，例如：平滑前：1111100000111100011平滑后：1111100000111111111解答：常规解法：[Copy to clipboard] [ - ]CODE:clear allI=[0 0 0 0 0 0 0 1 1 1 0 1 0 1 0 0 0 1 1 1 0 0 0 0 0 1 1 1 0 0 0 0 1 1 0 0 0 0 0 1];T = 4;Idiff = diff(I);index1 = find(Idiff == -1);index2 = find(Idiff == 1);if length(index1) == 0 | length(index2) == 0, %%%全1或者全0的情况，根据自己需要单独处理一下 coutinue;else if length(index1) == length(index2)+1, tt = index2 - index1(1:end-1); tt = [tt,length(I)-index1(end)]; index3 = find(tt <= T); for ii = 1:length(index3); I(index1(index3(ii))+1:index1(index3(ii))+tt(index3(ii))) = 1; end elseif length(index1) == length(index2)-1, if index2(1) <= T; I(1:index2(1)) = 1; end tt = [index2(2:end) - index1]; index3 = find(tt <= T); for ii = 1:length(index3); I(index1(index3(ii))+1:index1(index3(ii))+tt(index3(ii))) = 1; end elseif length(index1) == length(index2), if index1(1)>index2(1), if index2(1) <= T; I(1:index2(1)) = 1; end tt = index2(2:end) - index1(1:end-1); index3 = find(tt <= T); for ii = 1:length(index3); I(index1(index3(ii))+1:index1(index3(ii))+tt(index3(ii))) = 1; end if length(I)-index1(end) <=T, I(index1(end):end) = 1; end elseif index1(1) < index2(1), tt = index2 - index1; index3 = find(tt <= T); for ii = 1:length(index3); I(index1(index3(ii))+1:index1(index3(ii))+tt(index3(ii))) = 1; end end endendI优化解法：[Copy to clipboard] [ - ]CODE:a = [1 0 0 1 0 0 0 0 1 1 1 1 1 0 0 0 0 0 1 1 1 1 0 0 0 1 1 0 0 0 0 0 0 0 1 1 1 1 1 1];T = 4;b = sprintf('%d',a);[c,d]=regexp(b,'0+'); %%%ind = find(d - c < T);for ii = 1:length(ind), a(c(ind(ii)):d(ind(ii))) = 1;enda点评：优化解法巧妙地运用了正则表达式这个工具来处理特定字符串的搜索、匹配问题，因此比较快捷、简便。正则表达式的应用非常广泛，它是一种可以用于模式匹配和替换的强有力的工具。其作用如下：1) 测试字符串的某个模式。例如，可以对一个输入字符串进行测试，看在该字符串是否存在一个电话号码模式或一个信用卡号码模式。这称为数据有效性验证；2) 替换文本。可以在文档中使用一个正则表达式来标识特定文字，然后可以全部将其删除，或者替换为别的文字。备注：能否使用 regexpreg 命令一次性地进行搜索并完成替代呢？这样也许可以摆脱最后的循环语句。可惜本人功力有限，暂时未能实现。