计算机在生物工程中的应用

2MATLAB语言基础

2.1数据结构

MATLAB语言的赋值语句有两种：

变量名=运算表达式

[返回变量列表]=函数名(输入变量列表)

MATLAB支持变量和常量，更重要的，MATLAB支持IEEE标准的运算符

号，如pi为圆周率,i,j为虚数单位,Inf表示无穷大,NaN表示O/O、O*Inf或Inf/Inf

等运算结果，称为“非数”。MATLAB变量名应该由字母引导，后面可以跟数字、

字母或下划线等符号。MATLAB对变量名字母是区分大小写的。

2.1.1矩阵

MATLAB最基本的数据结构是复数矩阵。矩阵用方括号口输入，里面的各行

元素由分号分隔，而同一行的不同元素由逗号或空格分隔。如果赋值表达式末尾

有分号，则其结果将不显示，否则将显示出全部结果。

B=[l+9i,2+8i,3+7j;4+6j5+5i,6+4i;7+3i,8+2jli]

B=

1.0000+9.0000i2.0000+8.0000i3.0000+7.0000i

4.0000+6.0000i5.0000+5.0000i6.0000+4.0000i

7.0000+3.0000i8.0000+2.0000i0+l.OOOOi

其中，元素l+9i表示复数项。

MATLAB和其他语言不同，它无需事先声明矩阵的维数。下面的语句可以

建立一个更大的矩阵

B(2,5)=l

B=

1.0000+9.0000i2.0000+8.0000i3.0000+7.0000i00

4.0000+6.0000i5.0000+5.0000i6.0000+4.0000i01.0000

7.0000+3.0000i8.0000+2.0000i0+l.OOOOi00

(注意：B(2,5)代表第二行五列为1；与原来的对应，前面的儿个数不变，

新增加的后面几行其余的全部为零。)

对矩阵可以使用size。来测其大小，也可以使用reshape。函数重新按列排列。

对向量来说，还可以用length。来测其长度。不论原数组A是多少维的，A(:)将

返回列向量。

size(B)

length(B)%length函数代表“行”或“列”数中最大的一个。

reshape(B,5,3)%用此函数时元素个数不可改变，而且“重新排布”是按“列”

进行的。

B(：)

ans=

35

ans=

5

ans=

1.0000+9.0000i8.0000+2.0000i0

4.0000+6.0000i3.0000+7.0000i0

7.0000+3.0000i6.0000+4.0000i0

2.0000+8.0000i0+l.OOOOi1.0000

5.0000+5.0000i00

ans=

1.0000+9.0000i

4.0000+6.0000i

7.0000+3.0000i

2.0000+8.0000i

5.0000+5.0000i

8.0000+2.0000i

3.0000+7.0000i

6.0000+4.0000i

0+l.OOOOi

0

0

0

0

1.0000

0

2.1.2冒号表达式

冒号表达式：sl:s2:s3

冒号表达式是MATLAB里最具特色的表示方法。这一语句可以生成一个行

向量，其中si为向量的起始值，s2为步长，而s3为向量的终止值。例如

S=0:.l:2*pi

将产生一个起始于0,步距为0.1,而终止于6.2的向量，而不是终止于2*pi。如

果写成S=0:-0.l:2*pi;则不出现错误，而返回一-个空向量。

»s=(l:10)

12345678910

»s(l:-l:10)

ans=

Emptymatrix:l-by-0

另外，注意下列函数表达式及输出以作对比

»s(10:-l:l)

ans=

10987654321

»s(10:l:l)

ans=

Emptymatrix:1-by-0

冒号表达式可以用来提取矩阵元素，例如B(:,l)将提取B矩阵的第1列而

B(l:2,1:2:5)将提取B的前2行与1,3,5列组成的子矩阵。在矩阵提取时还可以采

用end这样的算符。如B(2:end,:)将提取B矩阵的后2行构成的子矩阵。

2.1.3其它数据结构

MATLAB的字符串是由单引号括起来的。如可以使用下面的命令赋值

strA=Thisisastring.1

MATLAB的单元数组可将复杂的数据结构纳入一个变量之下,单元数组由

大括号表示下标。

B={l；AlanShearer',180,[100,80,75;77,60,92;67,28,90;100,89,78]}

B=

[1]'AlanShearer'[180][4x3double]

访问单元数组应该用大括号进行，如第4单元中的元素可以由下面的语句得

出

B{4}

ans=

1008075

776092

672890

1008978

MATLAB的结构体有点象C语言的结构体数据结构。每个成员变量用点号

表示，如A.p表示A变量的p成员变量。用下面的语句可以建立一个小型的数

据库。

student_rec.number=l;

student_=,AlanShearer1;

student_rec.height=180;

student_rec.test=[100,80,75;77,60,92;67,28,90;100,89,78];

student_rec

student_rec=

number:1

name:'AlanShearer1

height:180

test:[4x3double]

2.2变量运算

2.2.1MATLAB变量的代数运算

如果给定两个矩阵A和B,则我们可以用A+B,A-B,A*B可以立即得出其加、

减和乘运算的结果。若这两个矩阵数学上不可以这样运算，则将得出错误信息，

并终止正在运行的程序。

在MATLAB下，如果A和B中有一个是标量，则可以无条件地进行这样的

运算。MATLAB不介意这些变量是纯实数还是含有虚部的复数。

矩阵的除法实际上就是线性方程的求解，如Ax=B这一线性方程组的解即

为x=inv(A)*B,inv(A)是A的逆矩阵。MATLAB中特有的左除\或右除/算法能更

有效地解这类线性方程组(包括矛盾方程和不定方程)，对矛盾方程会自动采用最

小二乘法求解。Ax=B的解用左除和右除分别表示为x=A\B或x'=(A\B)'=B7A',

号是矩阵转置符号。

K例2-21求解下面的四元一次方程组。

xl+x2+x3+x4=100

58.06x1+0.79x2+2.06x3+0.24x4=300

4.09x1+58.20x2+2.85x3+0.73x4=1171

20.77x1+29.99x2+43.26x3+4.91x4=1617

A=[1111

58.060.792.060.24

4.0958.202.850.73

20.7729.9943.264.91];

B=[100;300;1171;1617];

X=A\B

X=

4.1052

18.2863

15.7046

61.9039

方阵的乘方可以由“A”算符直接得出，如A^n。用MATLAB这样的语言，你

可以轻易地算出AA0.1,亦即A矩阵开10次方得出的主根。

矩阵的点运算也是相当重要的。所谓点运算即两个矩阵相应元素的运算，

如矩阵的点乘A.*B得出的是A和B对应元素的积，故一搬情况下A*B不等于

A.*B。点运算的概念又可以容易地用到点乘方上，例如A12,A.AA等都是可以接

受的运算式子。

2.2.2逻辑运算

MATLAB并没有单独定义逻辑变量。在MATLAB中，数值只有0和“非0”

的区分。非0往往被认为是逻辑真，或逻辑1。除了单独两个数值的逻辑运算外,

还支持矩阵的逻辑运算，如A&B,AIB＞hA分别表示逻辑与、或、非的运算。

请看下面的矩阵A和B与运算结果：

A=[0234;l350];B=[l053;1505];A&B

ans=

0011

1100

2.2.3关系表达式与表达式函数

MATLAB的大于、小于和等于等关系分别由〉、＜和==表示。判定方法不完

全等同于C语言这类只能处理单个标量的语言。MATLAB关系表达式返回的是

整个矩阵。例如，比较两个矩阵A和B是否相等，则可以给出如下命令，并得

出相应的结果

A=[0234;l35O];B=[1053;1505];A==B

ans=

0000

1000

确实使得A和B对应元素相等的位将返回1,否则返回0oMATLAB还可以用＞=

和＜=这样的符号来比较矩阵对应元素的大小。

另外，MATLAB还提供了all()和any()两个函数来对矩阵参数作逻辑判定。

all()函数在其中变元全部非0时返回1,而any()函数在变元有非零元素返回1。

find()函数将返回逻辑关系全部满足时的矩阵下标值，这个函数在编程中是相当

常用。

2.2.4其他运算

MATLAB还支持其他运算，如round。四舍五入取整、fix。向原点方向取整、

FLOOR。向数轴左边取整、CEILO向数轴右边取整、rem(X,Y)求X/Y的余数、

sign。返回符号等。

2.3流程控制

作为一种常用的编程语言，MATLAB支持各种流程控制结构，如循环结构、

条件转移结构、开关结构等另外MATLAB还支持一种新的结构…试探结构。

2.3.1循环结构

循环语句有两种结构：for...end结构和while…end结构。这两种语句结构

不完全相同，各有各的特色。for...end语句通常的调用格式为：

for循环变量=5193：52

循环体语句组

end

注意，这里的循环语句是以end结尾的，这和C语言的结构不完全一致。

K例2-32求1+2+...+100的值。

mysum=0;fori=1:1:100,mysum=mysum+i;end;mysum

mysum=

5050

在上面的式子中，可以看到for循环语句中S3的值为lo在MATLAB实际

编程中，如果S3的值为1,则可以在该语句中省略，故该语句可以简化成for

i=l:100o

在实际编程中，在MATLAB下采用循环语句会降低其执行速度，所以前面

的程序可以由内部求和函数sum()下面可提高执行速度：

i=l:l00;mysum=sum(i)或mysum=sum(l:100)

如果前面的100改成100000,再运行这一程序，则可以明显地看出，后一种

方法编写的程序比前一种方法快得多。测试执行某程序段所花费的时间可用命

令tic和toe实现，前者位于欲测试的程序段入口，后者位于出口。请看下面的

对比：

tic;mysum=0;fori=1:1:100000,mysum=mysum+i;end;mysum,toe

mysum=

5.000le+009

elapsed_time=

0.4110

tic;mysum=sum(1:100000),toe

mysum=

5.000le+009

elapsed_time=

0.0200

MATLAB并不要求循环点等间距，假设V为任意一个向量,则可以用fori=V

来表示循环：

V=[2.1-310/3piA2sin（pi/4）];mysum=0;fori=V,mysum=mysum+i;end;mysum

mysum=

13.0100

同样的问题在while循环结构下可以表示为

tic;mysum=0;i=1;while（i<=100000）,mysum=mysum+i;i=i+1;end;mysum,toc;

mysum=

5.000le+009

elapsed_time=

0.9110

由于while循环每次都要判断后面的条件式，并且循环内还多了i=i+l语句,

所以比for循环更慢！

2.3.2条件转移结构

条件转移语句和C语言相象

if条件式1

条件块语句组1

elseif条件式2

条件块语句组2

else（notes：此处else后不加IF表示如果前面没有满足则用

下面条件直接计算）

条件块语句组n+1

end

下面是利用条件转移语句把百分制转变为五级分制的例子：

mark=input。成绩='）;

ifmark>=90

disp（优）;

elseifmark>=80

disp。良）;

elseifmark>=70

disp（'中'）；

elseifmark>=60

dispC及格）;

else

dispC不及格）;

end

2.3.3开关结构

开关语句的基本结构为：

switch开关表达式

case表达式1

语句段1

case{表达式2,表达式3,...,表达式m}notes：只要满足括号中任意条件，均可执行下面的步骤

语句段2

otherwise

语句段n

end

MATLAB开关语句与C有区别：

无需像C语言那样在下一个case语句前加break语句；

把一些表达式用大括号括起来，中间用逗号分隔，当需满足表达式之一时执行某

一程序段；

各个表达式均不满足时，则将执行otherwise语句后面的语句段，此语句等价于

C语言中的default语句；

在case语句引导的各个表达式中，不要用重复的表达式，否则列在后面的开关

通路将永远也不能执行；

程序的执行结果和各个case语句的次序是无关的。

2.3.4试探结构

MATLAB从5.2版本开始提供了一种新的试探式语句结构，其一般的形式

为：

try

语句段1

catch

语句段2

end

本语句结构首先试探性地执行语句段1,如果在此段语句执行过程中出现错

误，则将错误信息赋给保留的lasterr变量，并放弃这段语句，转而执行语句段2

中的语句，这种新的语句结构是C等语言中所没有的。

try

a=sin(cos(pi/4);

catch

lasterr

end

???a=sin(cos(pi/4);

(此外输入至程序当中会显示错误的地方)

Error:expected,found.

2.4函数编写

MATLAB程序大致分为两类：M脚本文件(M-Script)和M函数(M-function),

它们均是普通的ASCII码构成的文件。M脚本文件中包含一族由MATLAB语

言所支持的语句，它类似于DOS下的批处理文件，它的执行方式很简单，用户

只需在MATLAB的提示符》下键入该M文件的文件名，这样MATLAB就会自

动执行该M文件中的各条语句，并将结果直接返回到MATLAB的工作空间。M

函数格式是MATLAB程序设计的主流，一般情况下，不建议您使用M脚本文件

格式编程。MATLAB提供了非常丰富的函数，利用菜单中的help命令可以了解

各种函数的功能和使用方法。

MATLAB的M函数是由function语句引导的，其基本格式如下：

function［返回变量列表］=函数名(输入变量列表)

注释说明语句段，由％引导-

输入、返回变量格式的检测

函数体语句

一(首词必须是function,注释那段作用为阐释愈义)

这里输入和返变量的实际个数分别由nargin和nargout两个MATLAB保

留变量来给出，只要进入该函数，MATLAB就将自动生成这两个变量，不论您

是否直接使用这两个变量。返回变量如果多于1个，则应该用方括号将它们括起

来，否则可以省去方括号。输入变量和返回变量之间用逗号来分割。

注释语句段的每行语句都应该由百分号引导，百分号后面的内容不执行，

只起注释作用。

一变量个数检测也是必要的，如果输入或返回变量格式不正确，则应该给出

相应的提示，不同的变量个数可能需要执行不同的程序段。

K例2.4』假设我们想生成一个nxm阶的Hilbert矩阵,它的第i行第j列的元素值

为l/(i+j-l)。我们想在编写的函数中实现下面儿点：

1.如果只给出一个输入参数，则会自动生成一个方阵，即令m=n

2.在函数中给出合适的帮助信息，包括基本功能、调用方式和参数说明

3.检测输入和返回变量的个数，如果有错误则给出错误信息

如果调用时不要求返回变量，则将显示结果矩阵。其实在编写程序时养成一

个好的习惯，无论对程序设计者还是对程序的维护者、使用者都是大有裨益的。

采用MATLAB函数编写格式和上述要求，我们可以编写出一个函数

functionA=myhilb(n,m)

%MYHILBademonstrativeM-function.

%A=MYHILB(N,M)generatesanNbyMHilbertmatrixA.

%A=MYHILB(N)generatesanNbyNsquareHilbertmatrix.

%MYHILB(N,M)displaysONLYtheHilbertmatrix,butdonotreturnany

%matrixbacktothecallingfunction.

%

%Seealso:HILB.

%DesignedbyProfessorDingyuXUE,NortheasternUniversity,PRC

%5April,1995,LastmodifiedbyDYXat21March,2000

ifnargout>1,error(Toomanyoutputarguments.1);end

ifnargin==l,m=n;

elseifnargin==0Inargin>2

error(Wrongnumberofiutputarguments.*);

end

Al=zeros(n,m);

fori=l:n

forj=l:m

Al(ij)=l/(i+j-l);

end,end

ifnargout==1,A=A1;elseifnargout==0,disp(Al);end

这样规范编写的函数用help命令可以显示出其帮助信息：

»helpmyhilb

MYHILBademonstrativeM-function.

A=MYHILB(N,M)generatesanNbyMHilbertmatrixA.

A=MYHILB(N)generatesanNbyNsquareHilbertmatrix.

MYHILB(N,M)displaysONLYtheHilbertmatrix,butdonotreturnany

matrixbacktothecallingfunction.

Seealso:HILB.

有了函数之后，可以采用下面的各种方法来调用它，并产生出所需的结果。

A=myhilb(3,4)

A=

1.00000.50000.33330.2500

0.50000.33330.25000.2000

0.33330.25000.20000.1667

A=myhilb(4)

A=

1.00000.50000.33330.2500

0.50000.33330.25000.2000

0.33330.25000.20000.1667

0.25000.20000.16670.1429

3科学计算的可视化

科学计算的可视化通过图形、图象等来表现，绘图是可视化核心内容。在传

统的编程语言中，绘图是编程中最繁杂的工作，而且费时费力。MATLAB提供

了非常方便实用的绘图功能，下面将介绍常用的绘图功能。

3.1二维绘图

科学计算的可视化通过图形、图象等来表现，绘图是可视化核心内容。在传

统的编程语言中，绘图是编程中最繁杂的工作，而且费时费力。MATLAB提供

了非常方便实用的绘图功能，下面将介绍常用的绘图功能。

假设有两个同长度的向量x和y,则用plot(x,y)就可以绘出二维曲线图。plot。

函数还可以同时绘制出多条曲线，如图3-l(a)(b)所示。

K例3-13绘制正弦余弦曲线:

t=0:.l:2*pi;

y=sin(t);

plot(t,y)

yl=cos(t);

figure

plot(t,y,t,yl);

如果未打开过图形窗口，将打开一个新的图形窗口绘图，若有打开过图形窗

口，则在最后使用的图形窗口上擦去原有图形并绘制新的图形。若想在另一个

新图形窗口中绘图，则在绘图前加一个figure语句，若想在一个指定的编号为n

的图形窗口中绘图，则在绘图前加一个figure(n)语句。

若两条曲线的横坐标向量相同，可以把两条曲线的纵坐标向量合并成一个矩

阵，每条曲线的纵坐标为矩阵的一行，绘图语句如下表示。

plot(t,[y;yl])

Figure2目回区

EileEditYiewInsertloolsDesktopWindowHelp

口自q昌ia|<CT)@|夏|口困|■回

（b）

图3-1二维曲线

axis。函数可以手动地设置x,y坐标轴范围。

axis（［02*pi-1.51.5］）

由MATLAB绘制的二维图形可以由下面的一些命令简单地修饰。如

grid%加网格线

xlabelC字符串）％给横坐标轴加说明

ylabelC字符串）％给纵坐标轴加说明，并自动旋转90度

title。字符串）％给整个图形加图题

plot。函数最完整的调用格式为：

pk>t（xl,yl,选项1,x2,y2,选项2,x3,y3,选项3,...）

其中所有的选项如表3-1所示。一些选项可连用，如才表示红色实线。

表3-1MATLAB绘图命令的各种选项

4.1:MATLAB绘图命令的各种选项

曲线线型曲线颜色标记符号

选项意义选项意义选项意义选项意义选项意

，一3实线,b,筮色'C'施绿色星号'pentagram't谆

>__>里1色，.，

虚线绿色'k'lii点号'o'网

卢线红紫色红色,x,叉T7号M,square,[

点划线5白色，y'黄色V'diacond'<

‘none'无线用一个1x3向俄任意指定，一，△,hexagram,六方

Lr,g,b］红球筮•:原色0<

用plotyy。函数可绘制具有两个纵坐标刻度的图形:

R例3-2』双纵坐标绘图：

t=0:.l:2*pi;

y=sin(t);

yl=sin(t).A2;

plotyy(t,y,t,yl);

图3-2双纵坐标绘图

3.2二维图形函数

除了标准的plot。函数外，MATLAB还提供了一些其他二维图形函数，具体调

用格式和意义请见下表

表3-2MATLAB提供的特殊二维绘图函数

表4.4:MATLAB提供的特殊二维曲线绘制函数

函数名意义常用调用格式

bar0.维条形图barCx,y)

comet()9星状轨迹图comet(x,y)

cocpass()罗薇图con-pass(x,y)

errorbar()误差限图形errorbar(x,y,1,u)

feather0羽毛状图feather(x,y)

fillO・维坤充函数fy,C)

hist()直方图hist(y,n)

loglog0对数图loglog(x,y)

polar()极坐标图polar(x,y)

quiver()磁力战图quiver(x,y)

stairs()阶梯图形stairs(x,y)

stem()火柒杆图stem(x9y)

semilogx()半对数图semilogx(x9y),semilogy(x,y)

极坐标曲线绘制：用polar(t,r)函数，其中t为角度向量，r为幅值向量。

K例3-33绘制t在［0,8］区间绘制r=cos(5t/4)+l/3的极坐标曲线。

t=0:.l:8*pi;

r=cos(5*t/4)+l/3;

polar(t,r)

图3-3极坐标绘图

坐标系的分割在MATLAB图形绘制中是很有特色的，分割方式的标准调用格式

为

subplot(n,m,k)

其中，n和m为将图形窗口分成的行数和列数，k为从左到右、从上到下循序的编

号。

K例3-4》在窗口上绘制四个图形：

x=-2:0.1:2;y=sin(x);

subplot(221);

feather(x,y);xlabel(*(a)feather。')

subplot(222);

stairs(x,y);xlabel(Xb)stairs()f)

subplot(223);

stem(x,y);xlabel(*(c)stem(),)

subplot(224);

fill(x,y,T);xlabel(Xd)fill。')

图3-4分区绘图

K例3-53考察MATLAB的Gauss伪随机数发生函数randn()的分布效果，先生成

30,000Gauss伪随机数，然后由hist()函数绘制出该伪随机数的分布函数，并和概率

密度的理论值

=L/x/SSe

相比较。

y=randn(1,30000);

xx=-3.8:0.4:3.8;

zz=hist(y,xx);

zz=zz/(30000*0.4);

bar(xx,zz)

xl=-3.8:0.1:3.8;

y1=l/sqrt(2*pi)*exp(-xl.A2/2);

holdon

plot(xl,yl)

holdoff

图3-5Gauss伪随机数概率分布图

holdon命令是为了原来图形窗口中补充绘图而不擦除原有的图形，绘完图后最

好马上用holdoff命令关闭。

可以使用半对数与全对数坐标系绘图，函数为：semilogx(),semilogyO和

loglogQo

K例3-63用半对数与全对数坐标系绘图：

theta=O:O.l:6*pi;r=cos(theta/3)+l/9;

subplot(2,2,l),polar(theta,r);

subplot(2,2,2);plot(theta,r);

subplot(2,2,3)；semilogx(theta,r);grid

subplot(2,2,4);semilogy(theta,r),grid

图3-6半对数坐标绘图

3.3三维图形函数

三维曲线绘图函数为plot3(),除了多了一个z坐标外，用法与二维曲线绘图函

数plot()类似。

K例3-71绘制x=sin(t),y=cos(t)和z=t螺旋线的一圈。

t=0:pi/50:2*pi;

x=sin(t);y=cos(t);z=t;

h=plot3(x,y,z,'g-');

set(h,'LineWidth',4*get(h,'LineWidth'));

gridon

图3-7三维曲线

三维网格图可以由mesh()函数绘制，其调用方法是mesh(x,y,z),其中x,y,z是

网格上的三坐标矩阵。通常x和y由meshgrid。函数生成。

K例3-81考虑下面给出的二元函数，在x,y平面内选择一个区域，然后绘制

出

z=/(工-2z)e'uJy

的三维表面图形。

首先可以调用meshgrid()函数生成x和y平面的网格表示。该函数的调用意义

十分明显，即可以产生一个横坐标起始于-3,中止于3,步距为0.1;纵坐标起始于

-2,中止于2,步距为0.1的网格分割。然后由公式计算出曲面的z矩阵。最后调

用mesh()函数绘制曲面的三维表面网格图形。

[x,y]=meshgrid(-3:0.1:3,-2:0.1:2);

z=(x.A2-2*x).*exp(-x.A2-y.A2-x.*y);

mesh(x,y,z)

图3-8三维表面图

MATLAB在三维绘图时，还可用其它绘图函数，如surf()、surfl()等，也有不同的

表面着色和插值方法。如shadingflat和shadinginterp两个命令将方便得出如下

的图形。

K例3-93surf()函数表面着色和插值举例。

surf(x,y,z),shadingflat

figure

surf(x,y,z),shadinginterp

图3-9网格与标注举例

3.4图形可视编辑

从MATLAB5.3版开始，提供了方便的图形编辑功能。在标准的MATLAB图

形窗口中有一个“图形编辑工具条“，典型窗口及编辑工具条如图3-10(a)所示。图

形窗口中提供了各种工具，允许用户自由地在图形上添加文字，箭头、线条等，

还允许用户任意地进行三维图的放大、缩小和旋转。

用户可以修改图形对象和坐标系的特征。如果想修改某个图形对象的特征，

则可以在编辑状态下，按下工具条中的选择(箭头)按钮，然后单击欲选对象，

被选中的对象上出现若干小黑方块的控制点时表示已被选中，双击对象可打开它

的属性编辑窗口，可对选中对象的属性进行修改，如图3-10所示。

图3-10编辑图形对象

单击空白处可选中坐标系，双击空白处可打开坐标系的属性编辑窗口，可对

坐标系的属性进行修改，如图3-ll(a)(b)所示。

图3-11编辑坐标系

用户还可以利用菜单中的功能进行编辑，如图3-12所示。

图3-12利用菜单功能编辑

3.5图象显示技术

351图象的读出与显示

图象的读出与显示要用到imread()和imshow()函数，程序如下:

I=imread('football.jpg');

imshow(I);

图3-13图象显示

3.5.2图象在柱面和上显示

图象在柱面上显示的关键函数是cylinder和warp()。举例如下:

I=imread(,football.jpg,);

[x,y,z]=cylinder;

warp(x,y,z,I);

图3-14图象在柱面上显示

3.5.3图象在球面上显示

图象在球面上显示的关键函数是sphere。和warp()0举例如下:

I=imread(,football.jpg,);

[x,y,z]=sphere(50);

warp(x,y,z,I);

图3-15图象在球面上显示

3.5.4图象在曲面上显示

图象在曲面上显示的关键首先计算出曲面的空间坐标，然后用warp()函数把

图象显示在曲面上即可。举例如下：

[x,y]=meshgrid(-3:0.1:3,-2:0.1:2);

z=(x.A2-2*x).*exp(-x.A2-y.A2-x.*y);

I=imread('football.jpg');

warp(x,y,z,I);

图3-16图象在曲面上显示

3.6动画与声音

凡是与时间有关的过程都是动态过程，动态过程可以用声像并茂的多媒体手

段来表现，如动画在多媒体教学过程中就得到广泛的应用。在矿物加工中，可以

通过监测磨机的声音来分析磨机的工作状况好坏，从而对磨机进行控制。

3.6.1动画的生成与播放

这里给出一个生成动画的程序，将要用到下面儿个关键函数peaks(),surfl(),

moviein(),view(),getframe。

I=imread('football.jpg');

[x,y,z]=peaks(30);

warp(x,y,z,I);

axis([-33-33-1010]);

axisoff;

n=80;

m=moviein(n);

fori=l:n

view(-42.5+4*(i-1),30+2*(i-1))

m(:,i)=getframe;

end

播放动画要用movie。函数。

movie(m)

图3-17动画演示

3.6.2声音的录制与播放

下面是声音的录制与播放程序,要用到wavrecord(),wavplay。函数。

fs=11025;

y=wavrecord(5*fs,fs/int16');

plot(y)

wavplay(y,fs);

8000

图3-18声音波形

画

囹4图象分析技术

图象分析在很多领域得到了越来越多的应用，在矿物加工领域也不例外。为

了让大家对图象分析有初步了解，下面给出图象分析中若干基本的分析。

关键函数：pixval,impixelinfo,imdistline,impixel,getline,improfile,imhist,

graythresh,im2bw,edge,imfill,imerode,imdilate,xor,line,frame2im,(imgetline,

imdivideline,improline,)bwarea,bweuler。

4.1位置尺寸测量

K例4-12用鼠标获取象素的位置与测量目标物体尺寸。

用pixval命令测量程序：

imshowrice.png;

pixval;

执行结果:

图4-1位置与尺寸测量

黑条中给出当前光标的坐标和灰度(对于彩色图象给出红绿蓝三个分量)，以及

鼠标拖出直线的长度。

与pixval命令功能相近的有测量位置颜色的impixelinfo函数和测量距离的

imdistline函数,举例如下：

imshowrice.png;

h=impixelinfo;

set(h；Position',[1503030020]);

imshowrice.png;

imdistline(gca,[121148],[110120]);

4.2选中象素坐标与颜色

K例4-2》显示选中象素的坐标和颜色。

用impixel函数测量程序:

imshowfootball.jpg;

[x,y,colors]=impixel

图4-2象素的颜色

执行结果：

x=

9

90

189

227

330

y=

83

85

83

82

79

colors=

246084

314460

245245243

15597109

NaNNaNNaN

注意:impixel的返回变量个数必须是0/或3,不能用2个返回变量,当只有一个返

回变量时,返回的是颜色而不是x坐标。若只需要返回坐标,还可用getline函数,

举例如下：

imshowfootball.jpg;

[x,y]=getline;

4.3折线象素坐标与颜色

K例4-31测量沿给定折线上像素的坐标与颜色。

improfile函数除了具有impixel和getline函数的功能外，还能返回折线上像素的坐

标和颜色，请看下面程序：

imshowfootball.jpg;

[CX,CY,C,x,y]=improfile

执行结果：

cx=

4

309

CY=

82

8

C(:,:,l)=

47

34

C(:,:,2)=

77

91

C(:,:,3)=

103

120

x=

4

91

91

309

y=

82

83

245

8

若没有返回变量,improfile将画出沿折线上的颜色分布图，如图4-3b所示。

improfile的返回变量可为0,1,3或5,当只有一个返回变量时，返回的是颜色。

可利用返回的直线坐标来切割图象，大家想想怎么做？

(a)给定折线

Figure2:Profile

(b)无返回变量时绘图输出

图4-3沿给定直线象素的灰度分布

4.4灰度直方图

K例4-41分析米粒的灰度图象，给出灰度直方图。

用imhist函数画灰度直方图程序：

I=imread('rice.png');

subplot(1,2,1);

imshow⑴;title('原图")

subplot(1,2,2);

imhist(I);title(，灰度直方图)

执行结果:

图4-4米粒的灰度直方图

灰度直方图右边的谷点是米粒与背景转变为黑白图象的最佳阈值，大约为140o

4.5转变为黑白图象

K例4-52把上题的米粒灰度图象转变为黑白及反白图象。

可用自动阈值函数graythresh函数和转换为黑白图象im2bw函数转换，负像只需

对正像取非运算即得，程序如下：

I=imread('rice.png');

level=graythresh(I);

%也可根据灰度直方图确定level,不用自动阈值函数graythresh

BW=im2bw(I,level);

subplot(1,2,1);

imshow(BW);titleC米粒的黑白图象)

BW=~im2bw(I,level);

subplot(1,2,2);

imshow(BW);title(，米粒的反白图象')

执行结果:

图4-5米粒的黑白及反白图象

4.6边缘检测

K例4-6》检测米粒灰度图象的米粒轮廓。

用edge函数检测米粒轮廓程序：

I=imread('rice.png');

BW=edge(I,'canny');

imshow(BW);

执行结果：

图4.6米粒轮廓图

edge函数的第二个输入参数为算法选择，有下面六种算法可选:

sobel

prewitt

roberts

log

zerocross

canny

不同的算法有不同的处理效果，算法的好坏与实际处理对象和目的有关，对于我

们处理的问题，canny是最好的。

4.7充填黑洞

K例4-73把上例米粒的图象边缘检测轮廓图中完整的米粒充填为白色。

用自动填补黑洞的imfil］函数继续处理上例的图象矩阵BW：

BWO=BW;BW=imfill(BW,'holes');imshow(BW);

执行结果：

图4-7自动填补黑洞后的图象

也可手工充填黑洞：

BW1=imfill(BWO);imshow(BW1);

%若充填不完全可重复执行上行继续充填

执行结果：

图4-8手工充填部分黑洞后的图象

4.8腐蚀与膨胀

在图4-7中一些不完整米粒的轮廓线尚未去掉，可以通过一次腐蚀imerode与

膨胀imdilate操作去掉。

这里有涉及到邻域的操作，我们只简单地介绍常用的4-邻域和8-邻域。在3X3

矩阵中心像素周围有1的位置都是中心像素的邻域。还可以自己定义其它的邻域

或更大的邻域，比如定义邻域为eye(5)0

4-邻域8-邻域

010111

111111

010111

腐蚀操作：

BW2=BW;

SE=[010;111;010];

BW=imerode(BW,SE);

imshow(BW);

执行结果:

图4-9腐蚀后的图象

膨胀操作：

BW=imdilate(BW,SE);

imshow(BW);

执行结果:

图4-10膨胀后的图象

经过一次腐蚀与膨胀操作，确实去掉了多余的轮廓线，但是要保留的完整米

粒是否能够在一次腐蚀和膨胀操作后恢复到原来的形状，可以用操作前后两图像

的异或xor操作来检查，结果表明用4-邻域是能够恢复，用8-邻域则恢复较差，

大家可以自行检验。

异或操作：

BW3=xor(BW,BW2);

imshow(BW3);

执行结果:

Figure1

图4-11膨胀后的图象

图4-10和图4-11实际上就是图4-7的分解为保留的部分和去掉的部分。

4.9图像分割

粘连在一起的米粒必须分割开，否则不能正确计数，可用前面介绍的improfile

函数获得折线上的坐标来实现分割，也可以用鼠标定位点坐标的来实现分割（两

点连线上各点的坐标需要插值），如improfile、impixel和getline函数都可用来

获得一系列的鼠标定位点。

若要保证在8-邻域下被认为是分割开的，分割线的像素总是沿上下左右方向

连续延伸（图4-12）,不能有跳跃。若有沿斜线方向延伸，在8-邻域下物体被认为

是联系在一起的，没有被真正分割开；但此时在4-邻域下被认为是分割开的。

100111111110111111

110011111111011111

111001111111101111

111100111111110111

111110011111111011

8-邻域割断4-邻域割断

Figure1EDI回区]

图4-12两种邻域下割断的分割线

8-邻域下，左图被识别为两个物体，右图在被识别为内部有3个黑洞的一个物体,

与边界相连的黑洞被认为是背景的部分，自动补洞只会充填不靠边界的三个小

黑洞，与边界相连的黑洞不会被充填。

设计程序要求可以连续进行分割，分割线的颜色取决于起点的颜色，当不需

要继续分割时，可在起点直接双击鼠标左键或单击鼠标右键退出(即只有一个坐

标点时退出)。

可用鼠标定位点的坐标来实现分割,getline、impixel和improfile都可以获得鼠

标定位点坐标。对于只保证在4-邻域下割断的场合，用line>frame2im和getframe

函数实现可以连续分割的简单程序如下。

BW4=BW;

imshow(BW);

[X,Y]=getline;

%[X,Y,C]=impixel;

%[CX,CY,C,X,Y]=improfile;

whilelength(X)>1

ifBW(round(Y(l)),round(X(l)))==0,C='k';else,C='w';end

line(X,Y,'Color',C)

MV=frame2im(getframe);

BW=MV(1:end-1,1:end-1,1)>0;

imshow(BW);

[X,Y]=getline;

end

如果原来的分割线旁偏移一个像素位置再画一条平行的分割线，就可以保证

在8-邻域下割断。下面自编的imgetline函数可通过第二个输入参数n实现4-邻

域或8-邻域割断的连续分割，n=4时为4-邻域割断，其它情况(包括缺省)皆为

8-邻域割断。

BW=imgetline(BW)

functionBW=imgetline(BW,n)

%IMGETLINEisafunctionofdrawinglinesinan%image.

%BW=IMGETLINE(BW,n)drawlinesinabinaryimage%BW,4-neibor

%regioncanbeselectedbyn=4,otherwise8-neiborregionisselected.

ifnargin-=2,n=8;elseifn〜=4,n=8;end

[X,Y]=getline;

whilelength(X)>1

ifBW(round(Y(l)),round(X(l)))==0,C=,k,;else,C=,w,;end

line(X,Y,Color',C)

ifn==8

Xl=X(2:end)-X(l:end-l);

Yl=Y(2:end)-Y(l:end-l);

YX=abs(Yl./Xl)>l;

fori=l:length(X)-l

ifYX(i)>l

Y(i)=Y(i)+l;

Y(i+l)=Y(i+l)+l;

else

X(i)=X(i)+l;

X(i+l)=X(i+l)+l;

end

end

line(X,Y；Color\C)

end

MV=frame2im(getframe);

BW=MV(1:end-1,1:end-1,1)>0;

imshow(BW);

fX,Y]=getline;

end

若在分割程序中不使用line>frame2im和getframe函数，也可以用下面自编

的imdivideline函数实现imgetline函数的功能，两者产生的分割线一般会有微小

差别。

BW=imdivideline(BW);

functionBW=imdivideline(BW,n)

%IMDIVIDELINEisafunctionofdrawinglinesinanimage.

%BW=IMDIVIDELINE(BW,n)drawlinesinabinaryimage%BW0,4-neibor

%regioncanbeselectedbyn=4,otherwise8-neiborregionisselected.

ifnargin-=2,n=8;elseifn〜=4,n=8;end

[X,Y]=getline;

whilelength(X)>l

L=l；

Xk=round(X(L));

Yk=round(Y(L));

C=BW(Yk,Xk);

whileL<length(X)

T=(Y(L+1)-Y(L))/(X(L+l)-X(L));

ifabs(T)>l

fork=round(Y(L)):sign(Y(L+l)-Y(L)):round(Y(L+l))

ifn==8,BW(k,Xk)=C;end

ifabs((k-Y(L))/T+X(L)-Xk)>0.5

Xk=Xk+sign(X(L+1)-Xk);

end

BW(k,Xk)=C;

end

else

fork=round(X(L)):sign(X(L+l)-X(L)):round(X(L+l))

ifn==8,BW(Yk,k)=C;end

ifabs((k-X(L))*T+Y(L)-Yk)>0.5

Yk=Yk+sign(Y(L+l)-Yk);

end

BW(Yk,k)=C;

end

end

L=L+1;

end

imshow(BW);

[X,Y]=getline;

end

下面是用improfile函数自编实现的4-邻域或8-邻域下割断的improline函数及

程序,第二个输入参数n=4时为4-邻域割断，其它情况(包括缺省)皆为8-邻

域割断，与前两个自编函数产生的分割线也会有微小差别。

BW=improline(BW);

functionBW=improline(BW,n)

%IMPROLINEisafunctionofdrawinglinesinan%image.

%BW=IMPROLINE(BW,n)drawlinesinabinaryimage%BW,4-neibor

%regioncanbeselectedbyn=4,otherwise8-neiborregionisselected.

ifnargin-=2,n=8;elseifn-=4,n=8;end

[CX,CY,C]=improfile;

whilelength(CX)>l

k=l;

Xk=round(CX(k));

Yk=round(CY(k));

C=BW(Yk,Xk);

whilek<length(CX)

k=k+1;

ifabs(CX(k)-Xk)>0.5

Xk=Xk+sign(CX(k)-Xk);

ifn==8,BW(Yk,Xk)=C;end

end

ifabs(CY(k)-Yk)>0.5

Yk=Yk+sign(CY(k)-Yk);

ifn==8,BW(Yk,Xk)=C;end

end

ifn==4,BW(Yk,Xk)=C;end

end

imshow(BW);

[CX,CY,C]=improfile;

end

执行结果:

图4-13分割粘连的米粒后的图象

4.10面积与计数

可以用bwarea函数统计黑白图象中白色物体占有的总像素数，用bweuler欧

拉数函数统计黑白图象中孤立白色物体的个数，它有第二个可选参数，为4时

按4-邻域割断识别独立的物体，为8时按8-邻域割断识别独立的物体，缺省值

为8o对于图4-12中4-邻域割断的图形,按8-邻域割断识别得到的欧拉数为-2

(一个物体减去三个黑洞)。由此就能计算出每个白色物体占有的平均像素数。

黑白图象的欧拉数定义为：

欧拉数=白色物体的个数-白色物体内黑色孔洞的个数

为了准确检测出白色物体的个数，事先必须填补白色物体内部的黑洞，用4-

邻域割断后不应再补洞，因为补洞会使分割线消失，而8-邻域割断后补洞不会影

响分割线。当白色物体内部没有黑洞时，欧拉数就等于白色物体个数。

程序：

%检测白色物体个数

N=bweuler(BW,4)

%检测白色区域总的象素个数

A=bwarea(BW)

%计算每个物体的平均象素个数

a=A/N

执行结果：

N=

70

A=

1.5447e+004

a=

220.6768

欧拉数函数的第一个输入参数为填补过黑洞的黑白图象矩阵，第二个参数为

邻域选择，可选为4和8对应于4-邻域和8-邻域,缺省为8-邻域。4-邻域比8-

邻域能更好区分粘连的颗粒。

[[例4-82检测米粒图像(文件名：rice.png)中米粒的个数及每个米粒的平均面积

(以像素为单位)。

下面程序是完全按照处理米粒的方法编写的程序，只是用了分区绘图并加了

文字说明，仍不够完善，大家可以修改得更好。

I=imread(Tice.png');

BW=edge(I,,canny,);

subplot(2,2,1);imshow(BW);

xlab&C边缘检测)

BW=imfill(BW,,holes,);

subplot(2,2,2);imshow(BW);

xlabelC自动填补黑洞)

BW2=BW;

SE=[010;l1l;010];

BW=imerode(BW,SE);

BW=imdilate(BW,SE);

subplot(2,2,3)；imshow(BW);

xlabelC腐蚀膨胀后图象)

BW3=xor(BW,BW2);

subplot(2,2,4);imshow(BW3);

xlabel(检查腐蚀膨胀效果)

figure;imshow(BW);

xlabelC分割粘连颗粒)

[X,Y]=getline;

whilelength(X)>l

ifBW(round(Y(1)),round(X(1)))==O,C='k';else,C='w';end

line(X,Y,'Color',C)

MV=frame2im(getframe);

BW=MV(l:end-l,l:end-1,1)>0;

imshow(BW);

fX,Y]=getline;

end

N=bweuler(BW,4)

A=bwarea(BW)

a=A/N

执行结果：

N=

70

A=

1.5450e+004

a=

220.7089

执行结果:

图4-14求米粒图像分析的中间过程

这是用米粒的处理程序处理的，这幅图像有两处需要分割，同时有一个米粒

轮廓线不封闭和两个与边界颗粒粘连的颗粒被去掉了，大家考虑能否修改程序，