計(jì)算機(jī)在生物工程中的應(yīng)用

2MATLAB語(yǔ)言基礎(chǔ)

2.1數(shù)據(jù)結(jié)構(gòu)

MATLAB語(yǔ)言的賦值語(yǔ)句有兩種：

變量名=運(yùn)算表達(dá)式

[返回變量列表]=函數(shù)名(輸入變量列表)

MATLAB支持變量和常量，更重要的，MATLAB支持IEEE標(biāo)準(zhǔn)的運(yùn)算符

號(hào)，如pi為圓周率,i,j為虛數(shù)單位,Inf表示無(wú)窮大,NaN表示O/O、O*Inf或Inf/Inf

等運(yùn)算結(jié)果，稱(chēng)為“非數(shù)”。MATLAB變量名應(yīng)該由字母引導(dǎo)，后面可以跟數(shù)字、

字母或下劃線等符號(hào)。MATLAB對(duì)變量名字母是區(qū)分大小寫(xiě)的。

2.1.1矩陣

MATLAB最基本的數(shù)據(jù)結(jié)構(gòu)是復(fù)數(shù)矩陣。矩陣用方括號(hào)口輸入，里面的各行

元素由分號(hào)分隔，而同一行的不同元素由逗號(hào)或空格分隔。如果賦值表達(dá)式末尾

有分號(hào)，則其結(jié)果將不顯示，否則將顯示出全部結(jié)果。

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表示復(fù)數(shù)項(xiàng)。

MATLAB和其他語(yǔ)言不同，它無(wú)需事先聲明矩陣的維數(shù)。下面的語(yǔ)句可以

建立一個(gè)更大的矩陣

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；與原來(lái)的對(duì)應(yīng)，前面的兒個(gè)數(shù)不變，

新增加的后面幾行其余的全部為零。)

對(duì)矩陣可以使用size。來(lái)測(cè)其大小，也可以使用reshape。函數(shù)重新按列排列。

對(duì)向量來(lái)說(shuō)，還可以用length。來(lái)測(cè)其長(zhǎng)度。不論原數(shù)組A是多少維的，A(:)將

返回列向量。

size(B)

length(B)%length函數(shù)代表“行”或“列”數(shù)中最大的一個(gè)。

reshape(B,5,3)%用此函數(shù)時(shí)元素個(gè)數(shù)不可改變，而且“重新排布”是按“列”

進(jìn)行的。

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冒號(hào)表達(dá)式

冒號(hào)表達(dá)式：sl:s2:s3

冒號(hào)表達(dá)式是MATLAB里最具特色的表示方法。這一語(yǔ)句可以生成一個(gè)行

向量，其中si為向量的起始值，s2為步長(zhǎng)，而s3為向量的終止值。例如

S=0:.l:2*pi

將產(chǎn)生一個(gè)起始于0,步距為0.1,而終止于6.2的向量，而不是終止于2*pi。如

果寫(xiě)成S=0:-0.l:2*pi;則不出現(xiàn)錯(cuò)誤，而返回一-個(gè)空向量。

?s=(l:10)

12345678910

?s(l:-l:10)

ans=

Emptymatrix:l-by-0

另外，注意下列函數(shù)表達(dá)式及輸出以作對(duì)比

?s(10:-l:l)

ans=

10987654321

?s(10:l:l)

ans=

Emptymatrix:1-by-0

冒號(hào)表達(dá)式可以用來(lái)提取矩陣元素，例如B(:,l)將提取B矩陣的第1列而

B(l:2,1:2:5)將提取B的前2行與1,3,5列組成的子矩陣。在矩陣提取時(shí)還可以采

用end這樣的算符。如B(2:end,:)將提取B矩陣的后2行構(gòu)成的子矩陣。

2.1.3其它數(shù)據(jù)結(jié)構(gòu)

MATLAB的字符串是由單引號(hào)括起來(lái)的。如可以使用下面的命令賦值

strA=Thisisastring.1

MATLAB的單元數(shù)組可將復(fù)雜的數(shù)據(jù)結(jié)構(gòu)納入一個(gè)變量之下,單元數(shù)組由

大括號(hào)表示下標(biāo)。

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

B=

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

訪問(wèn)單元數(shù)組應(yīng)該用大括號(hào)進(jìn)行，如第4單元中的元素可以由下面的語(yǔ)句得

出

B{4}

ans=

1008075

776092

672890

1008978

MATLAB的結(jié)構(gòu)體有點(diǎn)象C語(yǔ)言的結(jié)構(gòu)體數(shù)據(jù)結(jié)構(gòu)。每個(gè)成員變量用點(diǎn)號(hào)

表示，如A.p表示A變量的p成員變量。用下面的語(yǔ)句可以建立一個(gè)小型的數(shù)

據(jù)庫(kù)。

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變量運(yùn)算

2.2.1MATLAB變量的代數(shù)運(yùn)算

如果給定兩個(gè)矩陣A和B,則我們可以用A+B,A-B,A*B可以立即得出其加、

減和乘運(yùn)算的結(jié)果。若這兩個(gè)矩陣數(shù)學(xué)上不可以這樣運(yùn)算，則將得出錯(cuò)誤信息，

并終止正在運(yùn)行的程序。

在MATLAB下，如果A和B中有一個(gè)是標(biāo)量，則可以無(wú)條件地進(jìn)行這樣的

運(yùn)算。MATLAB不介意這些變量是純實(shí)數(shù)還是含有虛部的復(fù)數(shù)。

矩陣的除法實(shí)際上就是線性方程的求解，如Ax=B這一線性方程組的解即

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

有效地解這類(lèi)線性方程組(包括矛盾方程和不定方程)，對(duì)矛盾方程會(huì)自動(dòng)采用最

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

號(hào)是矩陣轉(zhuǎn)置符號(hào)。

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這樣的語(yǔ)言，你

可以輕易地算出AA0.1,亦即A矩陣開(kāi)10次方得出的主根。

矩陣的點(diǎn)運(yùn)算也是相當(dāng)重要的。所謂點(diǎn)運(yùn)算即兩個(gè)矩陣相應(yīng)元素的運(yùn)算，

如矩陣的點(diǎn)乘A.*B得出的是A和B對(duì)應(yīng)元素的積，故一搬情況下A*B不等于

A.*B。點(diǎn)運(yùn)算的概念又可以容易地用到點(diǎn)乘方上，例如A12,A.AA等都是可以接

受的運(yùn)算式子。

2.2.2邏輯運(yùn)算

MATLAB并沒(méi)有單獨(dú)定義邏輯變量。在MATLAB中，數(shù)值只有0和“非0”

的區(qū)分。非0往往被認(rèn)為是邏輯真，或邏輯1。除了單獨(dú)兩個(gè)數(shù)值的邏輯運(yùn)算外,

還支持矩陣的邏輯運(yùn)算，如A&B,AIB＞hA分別表示邏輯與、或、非的運(yùn)算。

請(qǐng)看下面的矩陣A和B與運(yùn)算結(jié)果：

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

ans=

0011

1100

2.2.3關(guān)系表達(dá)式與表達(dá)式函數(shù)

MATLAB的大于、小于和等于等關(guān)系分別由〉、＜和==表示。判定方法不完

全等同于C語(yǔ)言這類(lèi)只能處理單個(gè)標(biāo)量的語(yǔ)言。MATLAB關(guān)系表達(dá)式返回的是

整個(gè)矩陣。例如，比較兩個(gè)矩陣A和B是否相等，則可以給出如下命令，并得

出相應(yīng)的結(jié)果

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

ans=

0000

1000

確實(shí)使得A和B對(duì)應(yīng)元素相等的位將返回1,否則返回0oMATLAB還可以用＞=

和＜=這樣的符號(hào)來(lái)比較矩陣對(duì)應(yīng)元素的大小。

另外，MATLAB還提供了all()和any()兩個(gè)函數(shù)來(lái)對(duì)矩陣參數(shù)作邏輯判定。

all()函數(shù)在其中變?cè)糠?時(shí)返回1,而any()函數(shù)在變?cè)蟹橇阍胤祷?。

find()函數(shù)將返回邏輯關(guān)系全部滿(mǎn)足時(shí)的矩陣下標(biāo)值，這個(gè)函數(shù)在編程中是相當(dāng)

常用。

2.2.4其他運(yùn)算

MATLAB還支持其他運(yùn)算，如round。四舍五入取整、fix。向原點(diǎn)方向取整、

FLOOR。向數(shù)軸左邊取整、CEILO向數(shù)軸右邊取整、rem(X,Y)求X/Y的余數(shù)、

sign。返回符號(hào)等。

2.3流程控制

作為一種常用的編程語(yǔ)言，MATLAB支持各種流程控制結(jié)構(gòu)，如循環(huán)結(jié)構(gòu)、

條件轉(zhuǎn)移結(jié)構(gòu)、開(kāi)關(guān)結(jié)構(gòu)等另外MATLAB還支持一種新的結(jié)構(gòu)…試探結(jié)構(gòu)。

2.3.1循環(huán)結(jié)構(gòu)

循環(huán)語(yǔ)句有兩種結(jié)構(gòu)：for...end結(jié)構(gòu)和while…end結(jié)構(gòu)。這兩種語(yǔ)句結(jié)構(gòu)

不完全相同，各有各的特色。for...end語(yǔ)句通常的調(diào)用格式為：

for循環(huán)變量=5193：52

循環(huán)體語(yǔ)句組

end

注意，這里的循環(huán)語(yǔ)句是以end結(jié)尾的，這和C語(yǔ)言的結(jié)構(gòu)不完全一致。

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

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

mysum=

5050

在上面的式子中，可以看到for循環(huán)語(yǔ)句中S3的值為lo在MATLAB實(shí)際

編程中，如果S3的值為1,則可以在該語(yǔ)句中省略，故該語(yǔ)句可以簡(jiǎn)化成for

i=l:100o

在實(shí)際編程中，在MATLAB下采用循環(huán)語(yǔ)句會(huì)降低其執(zhí)行速度，所以前面

的程序可以由內(nèi)部求和函數(shù)sum()下面可提高執(zhí)行速度：

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

如果前面的100改成100000,再運(yùn)行這一程序，則可以明顯地看出，后一種

方法編寫(xiě)的程序比前一種方法快得多。測(cè)試執(zhí)行某程序段所花費(fèi)的時(shí)間可用命

令tic和toe實(shí)現(xiàn)，前者位于欲測(cè)試的程序段入口，后者位于出口。請(qǐng)看下面的

對(duì)比：

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并不要求循環(huán)點(diǎn)等間距，假設(shè)V為任意一個(gè)向量,則可以用fori=V

來(lái)表示循環(huán)：

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

mysum=

13.0100

同樣的問(wèn)題在while循環(huán)結(jié)構(gòu)下可以表示為

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循環(huán)每次都要判斷后面的條件式，并且循環(huán)內(nèi)還多了i=i+l語(yǔ)句,

所以比f(wàn)or循環(huán)更慢！

2.3.2條件轉(zhuǎn)移結(jié)構(gòu)

條件轉(zhuǎn)移語(yǔ)句和C語(yǔ)言相象

if條件式1

條件塊語(yǔ)句組1

elseif條件式2

條件塊語(yǔ)句組2

else（notes：此處else后不加IF表示如果前面沒(méi)有滿(mǎn)足則用

下面條件直接計(jì)算）

條件塊語(yǔ)句組n+1

end

下面是利用條件轉(zhuǎn)移語(yǔ)句把百分制轉(zhuǎn)變?yōu)槲寮?jí)分制的例子：

mark=input。成績(jī)='）;

ifmark>=90

disp（優(yōu)）;

elseifmark>=80

disp。良）;

elseifmark>=70

disp（'中'）；

elseifmark>=60

dispC及格）;

else

dispC不及格）;

end

2.3.3開(kāi)關(guān)結(jié)構(gòu)

開(kāi)關(guān)語(yǔ)句的基本結(jié)構(gòu)為：

switch開(kāi)關(guān)表達(dá)式

case表達(dá)式1

語(yǔ)句段1

case{表達(dá)式2,表達(dá)式3,...,表達(dá)式m}notes：只要滿(mǎn)足括號(hào)中任意條件，均可執(zhí)行下面的步驟

語(yǔ)句段2

otherwise

語(yǔ)句段n

end

MATLAB開(kāi)關(guān)語(yǔ)句與C有區(qū)別：

無(wú)需像C語(yǔ)言那樣在下一個(gè)case語(yǔ)句前加break語(yǔ)句；

把一些表達(dá)式用大括號(hào)括起來(lái)，中間用逗號(hào)分隔，當(dāng)需滿(mǎn)足表達(dá)式之一時(shí)執(zhí)行某

一程序段；

各個(gè)表達(dá)式均不滿(mǎn)足時(shí)，則將執(zhí)行otherwise語(yǔ)句后面的語(yǔ)句段，此語(yǔ)句等價(jià)于

C語(yǔ)言中的default語(yǔ)句；

在case語(yǔ)句引導(dǎo)的各個(gè)表達(dá)式中，不要用重復(fù)的表達(dá)式，否則列在后面的開(kāi)關(guān)

通路將永遠(yuǎn)也不能執(zhí)行；

程序的執(zhí)行結(jié)果和各個(gè)case語(yǔ)句的次序是無(wú)關(guān)的。

2.3.4試探結(jié)構(gòu)

MATLAB從5.2版本開(kāi)始提供了一種新的試探式語(yǔ)句結(jié)構(gòu)，其一般的形式

為：

try

語(yǔ)句段1

catch

語(yǔ)句段2

end

本語(yǔ)句結(jié)構(gòu)首先試探性地執(zhí)行語(yǔ)句段1,如果在此段語(yǔ)句執(zhí)行過(guò)程中出現(xiàn)錯(cuò)

誤，則將錯(cuò)誤信息賦給保留的lasterr變量，并放棄這段語(yǔ)句，轉(zhuǎn)而執(zhí)行語(yǔ)句段2

中的語(yǔ)句，這種新的語(yǔ)句結(jié)構(gòu)是C等語(yǔ)言中所沒(méi)有的。

try

a=sin(cos(pi/4);

catch

lasterr

end

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

(此外輸入至程序當(dāng)中會(huì)顯示錯(cuò)誤的地方)

Error:expected,found.

2.4函數(shù)編寫(xiě)

MATLAB程序大致分為兩類(lèi)：M腳本文件(M-Script)和M函數(shù)(M-function),

它們均是普通的ASCII碼構(gòu)成的文件。M腳本文件中包含一族由MATLAB語(yǔ)

言所支持的語(yǔ)句，它類(lèi)似于DOS下的批處理文件，它的執(zhí)行方式很簡(jiǎn)單，用戶(hù)

只需在MATLAB的提示符》下鍵入該M文件的文件名，這樣MATLAB就會(huì)自

動(dòng)執(zhí)行該M文件中的各條語(yǔ)句，并將結(jié)果直接返回到MATLAB的工作空間。M

函數(shù)格式是MATLAB程序設(shè)計(jì)的主流，一般情況下，不建議您使用M腳本文件

格式編程。MATLAB提供了非常豐富的函數(shù)，利用菜單中的help命令可以了解

各種函數(shù)的功能和使用方法。

MATLAB的M函數(shù)是由function語(yǔ)句引導(dǎo)的，其基本格式如下：

function［返回變量列表］=函數(shù)名(輸入變量列表)

注釋說(shuō)明語(yǔ)句段，由％引導(dǎo)-

輸入、返回變量格式的檢測(cè)

函數(shù)體語(yǔ)句

一(首詞必須是function,注釋那段作用為闡釋愈義)

這里輸入和返變量的實(shí)際個(gè)數(shù)分別由nargin和nargout兩個(gè)MATLAB保

留變量來(lái)給出，只要進(jìn)入該函數(shù)，MATLAB就將自動(dòng)生成這兩個(gè)變量，不論您

是否直接使用這兩個(gè)變量。返回變量如果多于1個(gè)，則應(yīng)該用方括號(hào)將它們括起

來(lái)，否則可以省去方括號(hào)。輸入變量和返回變量之間用逗號(hào)來(lái)分割。

注釋語(yǔ)句段的每行語(yǔ)句都應(yīng)該由百分號(hào)引導(dǎo)，百分號(hào)后面的內(nèi)容不執(zhí)行，

只起注釋作用。

一變量個(gè)數(shù)檢測(cè)也是必要的，如果輸入或返回變量格式不正確，則應(yīng)該給出

相應(yīng)的提示，不同的變量個(gè)數(shù)可能需要執(zhí)行不同的程序段。

K例2.4』假設(shè)我們想生成一個(gè)nxm階的Hilbert矩陣,它的第i行第j列的元素值

為l/(i+j-l)。我們想在編寫(xiě)的函數(shù)中實(shí)現(xiàn)下面兒點(diǎn)：

1.如果只給出一個(gè)輸入?yún)?shù)，則會(huì)自動(dòng)生成一個(gè)方陣，即令m=n

2.在函數(shù)中給出合適的幫助信息，包括基本功能、調(diào)用方式和參數(shù)說(shuō)明

3.檢測(cè)輸入和返回變量的個(gè)數(shù)，如果有錯(cuò)誤則給出錯(cuò)誤信息

如果調(diào)用時(shí)不要求返回變量，則將顯示結(jié)果矩陣。其實(shí)在編寫(xiě)程序時(shí)養(yǎng)成一

個(gè)好的習(xí)慣，無(wú)論對(duì)程序設(shè)計(jì)者還是對(duì)程序的維護(hù)者、使用者都是大有裨益的。

采用MATLAB函數(shù)編寫(xiě)格式和上述要求，我們可以編寫(xiě)出一個(gè)函數(shù)

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

這樣規(guī)范編寫(xiě)的函數(shù)用help命令可以顯示出其幫助信息：

?helpmyhilb

MYHILBademonstrativeM-function.

A=MYHILB(N,M)generatesanNbyMHilbertmatrixA.

A=MYHILB(N)generatesanNbyNsquareHilbertmatrix.

MYHILB(N,M)displaysONLYtheHilbertmatrix,butdonotreturnany

matrixbacktothecallingfunction.

Seealso:HILB.

有了函數(shù)之后，可以采用下面的各種方法來(lái)調(diào)用它，并產(chǎn)生出所需的結(jié)果。

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科學(xué)計(jì)算的可視化

科學(xué)計(jì)算的可視化通過(guò)圖形、圖象等來(lái)表現(xiàn)，繪圖是可視化核心內(nèi)容。在傳

統(tǒng)的編程語(yǔ)言中，繪圖是編程中最繁雜的工作，而且費(fèi)時(shí)費(fèi)力。MATLAB提供

了非常方便實(shí)用的繪圖功能，下面將介紹常用的繪圖功能。

3.1二維繪圖

科學(xué)計(jì)算的可視化通過(guò)圖形、圖象等來(lái)表現(xiàn)，繪圖是可視化核心內(nèi)容。在傳

統(tǒng)的編程語(yǔ)言中，繪圖是編程中最繁雜的工作，而且費(fèi)時(shí)費(fèi)力。MATLAB提供

了非常方便實(shí)用的繪圖功能，下面將介紹常用的繪圖功能。

假設(shè)有兩個(gè)同長(zhǎng)度的向量x和y,則用plot(x,y)就可以繪出二維曲線圖。plot。

函數(shù)還可以同時(shí)繪制出多條曲線，如圖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);

如果未打開(kāi)過(guò)圖形窗口，將打開(kāi)一個(gè)新的圖形窗口繪圖，若有打開(kāi)過(guò)圖形窗

口，則在最后使用的圖形窗口上擦去原有圖形并繪制新的圖形。若想在另一個(gè)

新圖形窗口中繪圖，則在繪圖前加一個(gè)figure語(yǔ)句，若想在一個(gè)指定的編號(hào)為n

的圖形窗口中繪圖，則在繪圖前加一個(gè)figure(n)語(yǔ)句。

若兩條曲線的橫坐標(biāo)向量相同，可以把兩條曲線的縱坐標(biāo)向量合并成一個(gè)矩

陣，每條曲線的縱坐標(biāo)為矩陣的一行，繪圖語(yǔ)句如下表示。

plot(t,[y;yl])

Figure2目回區(qū)

EileEditYiewInsertloolsDesktopWindowHelp

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

（b）

圖3-1二維曲線

axis。函數(shù)可以手動(dòng)地設(shè)置x,y坐標(biāo)軸范圍。

axis（［02*pi-1.51.5］）

由MATLAB繪制的二維圖形可以由下面的一些命令簡(jiǎn)單地修飾。如

grid%加網(wǎng)格線

xlabelC字符串）％給橫坐標(biāo)軸加說(shuō)明

ylabelC字符串）％給縱坐標(biāo)軸加說(shuō)明，并自動(dòng)旋轉(zhuǎn)90度

title。字符串）％給整個(gè)圖形加圖題

plot。函數(shù)最完整的調(diào)用格式為：

pk>t（xl,yl,選項(xiàng)1,x2,y2,選項(xiàng)2,x3,y3,選項(xiàng)3,...）

其中所有的選項(xiàng)如表3-1所示。一些選項(xiàng)可連用，如才表示紅色實(shí)線。

表3-1MATLAB繪圖命令的各種選項(xiàng)

4.1:MATLAB繪圖命令的各種選項(xiàng)

曲線線型曲線顏色標(biāo)記符號(hào)

選項(xiàng)意義選項(xiàng)意義選項(xiàng)意義選項(xiàng)意義選項(xiàng)意

，一3實(shí)線,b,筮色'C'施綠色星號(hào)'pentagram't諄

>__>里1色，.，

虛線綠色'k'lii點(diǎn)號(hào)'o'網(wǎng)

盧線紅紫色紅色,x,叉T7號(hào)M,square,[

點(diǎn)劃線5白色，y'黃色V'diacond'<

‘none'無(wú)線用一個(gè)1x3向俄任意指定，一，△,hexagram,六方

Lr,g,b］紅球筮?:原色0<

用plotyy。函數(shù)可繪制具有兩個(gè)縱坐標(biāo)刻度的圖形:

R例3-2』雙縱坐標(biāo)繪圖：

t=0:.l:2*pi;

y=sin(t);

yl=sin(t).A2;

plotyy(t,y,t,yl);

圖3-2雙縱坐標(biāo)繪圖

3.2二維圖形函數(shù)

除了標(biāo)準(zhǔn)的plot。函數(shù)外，MATLAB還提供了一些其他二維圖形函數(shù)，具體調(diào)

用格式和意義請(qǐng)見(jiàn)下表

表3-2MATLAB提供的特殊二維繪圖函數(shù)

表4.4:MATLAB提供的特殊二維曲線繪制函數(shù)

函數(shù)名意義常用調(diào)用格式

bar0.維條形圖barCx,y)

comet()9星狀軌跡圖comet(x,y)

cocpass()羅薇圖con-pass(x,y)

errorbar()誤差限圖形errorbar(x,y,1,u)

feather0羽毛狀圖feather(x,y)

fillO?維坤充函數(shù)fy,C)

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

loglog0對(duì)數(shù)圖loglog(x,y)

polar()極坐標(biāo)圖polar(x,y)

quiver()磁力戰(zhàn)圖quiver(x,y)

stairs()階梯圖形stairs(x,y)

stem()火柒桿圖stem(x9y)

semilogx()半對(duì)數(shù)圖semilogx(x9y),semilogy(x,y)

極坐標(biāo)曲線繪制：用polar(t,r)函數(shù)，其中t為角度向量，r為幅值向量。

K例3-33繪制t在［0,8］區(qū)間繪制r=cos(5t/4)+l/3的極坐標(biāo)曲線。

t=0:.l:8*pi;

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

polar(t,r)

圖3-3極坐標(biāo)繪圖

坐標(biāo)系的分割在MATLAB圖形繪制中是很有特色的，分割方式的標(biāo)準(zhǔn)調(diào)用格式

為

subplot(n,m,k)

其中，n和m為將圖形窗口分成的行數(shù)和列數(shù)，k為從左到右、從上到下循序的編

號(hào)。

K例3-4》在窗口上繪制四個(gè)圖形：

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分區(qū)繪圖

K例3-53考察MATLAB的Gauss偽隨機(jī)數(shù)發(fā)生函數(shù)randn()的分布效果，先生成

30,000Gauss偽隨機(jī)數(shù)，然后由hist()函數(shù)繪制出該偽隨機(jī)數(shù)的分布函數(shù)，并和概率

密度的理論值

=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偽隨機(jī)數(shù)概率分布圖

holdon命令是為了原來(lái)圖形窗口中補(bǔ)充繪圖而不擦除原有的圖形，繪完圖后最

好馬上用holdoff命令關(guān)閉。

可以使用半對(duì)數(shù)與全對(duì)數(shù)坐標(biāo)系繪圖，函數(shù)為：semilogx(),semilogyO和

loglogQo

K例3-63用半對(duì)數(shù)與全對(duì)數(shù)坐標(biāo)系繪圖：

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半對(duì)數(shù)坐標(biāo)繪圖

3.3三維圖形函數(shù)

三維曲線繪圖函數(shù)為plot3(),除了多了一個(gè)z坐標(biāo)外，用法與二維曲線繪圖函

數(shù)plot()類(lèi)似。

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三維曲線

三維網(wǎng)格圖可以由mesh()函數(shù)繪制，其調(diào)用方法是mesh(x,y,z),其中x,y,z是

網(wǎng)格上的三坐標(biāo)矩陣。通常x和y由meshgrid。函數(shù)生成。

K例3-81考慮下面給出的二元函數(shù)，在x,y平面內(nèi)選擇一個(gè)區(qū)域，然后繪制

出

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

的三維表面圖形。

首先可以調(diào)用meshgrid()函數(shù)生成x和y平面的網(wǎng)格表示。該函數(shù)的調(diào)用意義

十分明顯，即可以產(chǎn)生一個(gè)橫坐標(biāo)起始于-3,中止于3,步距為0.1;縱坐標(biāo)起始于

-2,中止于2,步距為0.1的網(wǎng)格分割。然后由公式計(jì)算出曲面的z矩陣。最后調(diào)

用mesh()函數(shù)繪制曲面的三維表面網(wǎng)格圖形。

[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在三維繪圖時(shí)，還可用其它繪圖函數(shù)，如surf()、surfl()等，也有不同的

表面著色和插值方法。如shadingflat和shadinginterp兩個(gè)命令將方便得出如下

的圖形。

K例3-93surf()函數(shù)表面著色和插值舉例。

surf(x,y,z),shadingflat

figure

surf(x,y,z),shadinginterp

圖3-9網(wǎng)格與標(biāo)注舉例

3.4圖形可視編輯

從MATLAB5.3版開(kāi)始，提供了方便的圖形編輯功能。在標(biāo)準(zhǔn)的MATLAB圖

形窗口中有一個(gè)“圖形編輯工具條“，典型窗口及編輯工具條如圖3-10(a)所示。圖

形窗口中提供了各種工具，允許用戶(hù)自由地在圖形上添加文字，箭頭、線條等，

還允許用戶(hù)任意地進(jìn)行三維圖的放大、縮小和旋轉(zhuǎn)。

用戶(hù)可以修改圖形對(duì)象和坐標(biāo)系的特征。如果想修改某個(gè)圖形對(duì)象的特征，

則可以在編輯狀態(tài)下，按下工具條中的選擇(箭頭)按鈕，然后單擊欲選對(duì)象，

被選中的對(duì)象上出現(xiàn)若干小黑方塊的控制點(diǎn)時(shí)表示已被選中，雙擊對(duì)象可打開(kāi)它

的屬性編輯窗口，可對(duì)選中對(duì)象的屬性進(jìn)行修改，如圖3-10所示。

圖3-10編輯圖形對(duì)象

單擊空白處可選中坐標(biāo)系，雙擊空白處可打開(kāi)坐標(biāo)系的屬性編輯窗口，可對(duì)

坐標(biāo)系的屬性進(jìn)行修改，如圖3-ll(a)(b)所示。

圖3-11編輯坐標(biāo)系

用戶(hù)還可以利用菜單中的功能進(jìn)行編輯，如圖3-12所示。

圖3-12利用菜單功能編輯

3.5圖象顯示技術(shù)

351圖象的讀出與顯示

圖象的讀出與顯示要用到imread()和imshow()函數(shù)，程序如下:

I=imread('football.jpg');

imshow(I);

圖3-13圖象顯示

3.5.2圖象在柱面和上顯示

圖象在柱面上顯示的關(guān)鍵函數(shù)是cylinder和warp()。舉例如下:

I=imread(,football.jpg,);

[x,y,z]=cylinder;

warp(x,y,z,I);

圖3-14圖象在柱面上顯示

3.5.3圖象在球面上顯示

圖象在球面上顯示的關(guān)鍵函數(shù)是sphere。和warp()0舉例如下:

I=imread(,football.jpg,);

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

warp(x,y,z,I);

圖3-15圖象在球面上顯示

3.5.4圖象在曲面上顯示

圖象在曲面上顯示的關(guān)鍵首先計(jì)算出曲面的空間坐標(biāo)，然后用warp()函數(shù)把

圖象顯示在曲面上即可。舉例如下：

[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動(dòng)畫(huà)與聲音

凡是與時(shí)間有關(guān)的過(guò)程都是動(dòng)態(tài)過(guò)程，動(dòng)態(tài)過(guò)程可以用聲像并茂的多媒體手

段來(lái)表現(xiàn)，如動(dòng)畫(huà)在多媒體教學(xué)過(guò)程中就得到廣泛的應(yīng)用。在礦物加工中，可以

通過(guò)監(jiān)測(cè)磨機(jī)的聲音來(lái)分析磨機(jī)的工作狀況好壞，從而對(duì)磨機(jī)進(jìn)行控制。

3.6.1動(dòng)畫(huà)的生成與播放

這里給出一個(gè)生成動(dòng)畫(huà)的程序，將要用到下面兒個(gè)關(guān)鍵函數(shù)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

播放動(dòng)畫(huà)要用movie。函數(shù)。

movie(m)

圖3-17動(dòng)畫(huà)演示

3.6.2聲音的錄制與播放

下面是聲音的錄制與播放程序,要用到wavrecord(),wavplay。函數(shù)。

fs=11025;

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

plot(y)

wavplay(y,fs);

8000

圖3-18聲音波形

畫(huà)

囹4圖象分析技術(shù)

圖象分析在很多領(lǐng)域得到了越來(lái)越多的應(yīng)用，在礦物加工領(lǐng)域也不例外。為

了讓大家對(duì)圖象分析有初步了解，下面給出圖象分析中若干基本的分析。

關(guān)鍵函數(shù)：pixval,impixelinfo,imdistline,impixel,getline,improfile,imhist,

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

imdivideline,improline,)bwarea,bweuler。

4.1位置尺寸測(cè)量

K例4-12用鼠標(biāo)獲取象素的位置與測(cè)量目標(biāo)物體尺寸。

用pixval命令測(cè)量程序：

imshowrice.png;

pixval;

執(zhí)行結(jié)果:

圖4-1位置與尺寸測(cè)量

黑條中給出當(dāng)前光標(biāo)的坐標(biāo)和灰度(對(duì)于彩色圖象給出紅綠藍(lán)三個(gè)分量)，以及

鼠標(biāo)拖出直線的長(zhǎng)度。

與pixval命令功能相近的有測(cè)量位置顏色的impixelinfo函數(shù)和測(cè)量距離的

imdistline函數(shù),舉例如下：

imshowrice.png;

h=impixelinfo;

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

imshowrice.png;

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

4.2選中象素坐標(biāo)與顏色

K例4-2》顯示選中象素的坐標(biāo)和顏色。

用impixel函數(shù)測(cè)量程序:

imshowfootball.jpg;

[x,y,colors]=impixel

圖4-2象素的顏色

執(zhí)行結(jié)果：

x=

9

90

189

227

330

y=

83

85

83

82

79

colors=

246084

314460

245245243

15597109

NaNNaNNaN

注意:impixel的返回變量個(gè)數(shù)必須是0/或3,不能用2個(gè)返回變量,當(dāng)只有一個(gè)返

回變量時(shí),返回的是顏色而不是x坐標(biāo)。若只需要返回坐標(biāo),還可用getline函數(shù),

舉例如下：

imshowfootball.jpg;

[x,y]=getline;

4.3折線象素坐標(biāo)與顏色

K例4-31測(cè)量沿給定折線上像素的坐標(biāo)與顏色。

improfile函數(shù)除了具有impixel和getline函數(shù)的功能外，還能返回折線上像素的坐

標(biāo)和顏色，請(qǐng)看下面程序：

imshowfootball.jpg;

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

執(zhí)行結(jié)果：

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

若沒(méi)有返回變量,improfile將畫(huà)出沿折線上的顏色分布圖，如圖4-3b所示。

improfile的返回變量可為0,1,3或5,當(dāng)只有一個(gè)返回變量時(shí)，返回的是顏色。

可利用返回的直線坐標(biāo)來(lái)切割圖象，大家想想怎么做？

(a)給定折線

Figure2:Profile

(b)無(wú)返回變量時(shí)繪圖輸出

圖4-3沿給定直線象素的灰度分布

4.4灰度直方圖

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

用imhist函數(shù)畫(huà)灰度直方圖程序：

I=imread('rice.png');

subplot(1,2,1);

imshow⑴;title('原圖")

subplot(1,2,2);

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

執(zhí)行結(jié)果:

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

灰度直方圖右邊的谷點(diǎn)是米粒與背景轉(zhuǎn)變?yōu)楹诎讏D象的最佳閾值，大約為140o

4.5轉(zhuǎn)變?yōu)楹诎讏D象

K例4-52把上題的米粒灰度圖象轉(zhuǎn)變?yōu)楹诎准胺窗讏D象。

可用自動(dòng)閾值函數(shù)graythresh函數(shù)和轉(zhuǎn)換為黑白圖象im2bw函數(shù)轉(zhuǎn)換，負(fù)像只需

對(duì)正像取非運(yùn)算即得，程序如下：

I=imread('rice.png');

level=graythresh(I);

%也可根據(jù)灰度直方圖確定level,不用自動(dòng)閾值函數(shù)graythresh

BW=im2bw(I,level);

subplot(1,2,1);

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

BW=~im2bw(I,level);

subplot(1,2,2);

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

執(zhí)行結(jié)果:

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

4.6邊緣檢測(cè)

K例4-6》檢測(cè)米?；叶葓D象的米粒輪廓。

用edge函數(shù)檢測(cè)米粒輪廓程序：

I=imread('rice.png');

BW=edge(I,'canny');

imshow(BW);

執(zhí)行結(jié)果：

圖4.6米粒輪廓圖

edge函數(shù)的第二個(gè)輸入?yún)?shù)為算法選擇，有下面六種算法可選:

sobel

prewitt

roberts

log

zerocross

canny

不同的算法有不同的處理效果，算法的好壞與實(shí)際處理對(duì)象和目的有關(guān)，對(duì)于我

們處理的問(wèn)題，canny是最好的。

4.7充填黑洞

K例4-73把上例米粒的圖象邊緣檢測(cè)輪廓圖中完整的米粒充填為白色。

用自動(dòng)填補(bǔ)黑洞的imfil］函數(shù)繼續(xù)處理上例的圖象矩陣BW：

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

執(zhí)行結(jié)果：

圖4-7自動(dòng)填補(bǔ)黑洞后的圖象

也可手工充填黑洞：

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

%若充填不完全可重復(fù)執(zhí)行上行繼續(xù)充填

執(zhí)行結(jié)果：

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

4.8腐蝕與膨脹

在圖4-7中一些不完整米粒的輪廓線尚未去掉，可以通過(guò)一次腐蝕imerode與

膨脹imdilate操作去掉。

這里有涉及到鄰域的操作，我們只簡(jiǎn)單地介紹常用的4-鄰域和8-鄰域。在3X3

矩陣中心像素周?chē)?的位置都是中心像素的鄰域。還可以自己定義其它的鄰域

或更大的鄰域，比如定義鄰域?yàn)閑ye(5)0

4-鄰域8-鄰域

010111

111111

010111

腐蝕操作：

BW2=BW;

SE=[010;111;010];

BW=imerode(BW,SE);

imshow(BW);

執(zhí)行結(jié)果:

圖4-9腐蝕后的圖象

膨脹操作：

BW=imdilate(BW,SE);

imshow(BW);

執(zhí)行結(jié)果:

圖4-10膨脹后的圖象

經(jīng)過(guò)一次腐蝕與膨脹操作，確實(shí)去掉了多余的輪廓線，但是要保留的完整米

粒是否能夠在一次腐蝕和膨脹操作后恢復(fù)到原來(lái)的形狀，可以用操作前后兩圖像

的異或xor操作來(lái)檢查，結(jié)果表明用4-鄰域是能夠恢復(fù)，用8-鄰域則恢復(fù)較差，

大家可以自行檢驗(yàn)。

異或操作：

BW3=xor(BW,BW2);

imshow(BW3);

執(zhí)行結(jié)果:

Figure1

圖4-11膨脹后的圖象

圖4-10和圖4-11實(shí)際上就是圖4-7的分解為保留的部分和去掉的部分。

4.9圖像分割

粘連在一起的米粒必須分割開(kāi)，否則不能正確計(jì)數(shù)，可用前面介紹的improfile

函數(shù)獲得折線上的坐標(biāo)來(lái)實(shí)現(xiàn)分割，也可以用鼠標(biāo)定位點(diǎn)坐標(biāo)的來(lái)實(shí)現(xiàn)分割（兩

點(diǎn)連線上各點(diǎn)的坐標(biāo)需要插值），如improfile、impixel和getline函數(shù)都可用來(lái)

獲得一系列的鼠標(biāo)定位點(diǎn)。

若要保證在8-鄰域下被認(rèn)為是分割開(kāi)的，分割線的像素總是沿上下左右方向

連續(xù)延伸（圖4-12）,不能有跳躍。若有沿斜線方向延伸，在8-鄰域下物體被認(rèn)為

是聯(lián)系在一起的，沒(méi)有被真正分割開(kāi)；但此時(shí)在4-鄰域下被認(rèn)為是分割開(kāi)的。

100111111110111111

110011111111011111

111001111111101111

111100111111110111

111110011111111011

8-鄰域割斷4-鄰域割斷

Figure1EDI回區(qū)]

圖4-12兩種鄰域下割斷的分割線

8-鄰域下，左圖被識(shí)別為兩個(gè)物體，右圖在被識(shí)別為內(nèi)部有3個(gè)黑洞的一個(gè)物體,

與邊界相連的黑洞被認(rèn)為是背景的部分，自動(dòng)補(bǔ)洞只會(huì)充填不靠邊界的三個(gè)小

黑洞，與邊界相連的黑洞不會(huì)被充填。

設(shè)計(jì)程序要求可以連續(xù)進(jìn)行分割，分割線的顏色取決于起點(diǎn)的顏色，當(dāng)不需

要繼續(xù)分割時(shí)，可在起點(diǎn)直接雙擊鼠標(biāo)左鍵或單擊鼠標(biāo)右鍵退出(即只有一個(gè)坐

標(biāo)點(diǎn)時(shí)退出)。

可用鼠標(biāo)定位點(diǎn)的坐標(biāo)來(lái)實(shí)現(xiàn)分割,getline、impixel和improfile都可以獲得鼠

標(biāo)定位點(diǎn)坐標(biāo)。對(duì)于只保證在4-鄰域下割斷的場(chǎng)合，用line>frame2im和getframe

函數(shù)實(shí)現(xiàn)可以連續(xù)分割的簡(jiǎn)單程序如下。

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

如果原來(lái)的分割線旁偏移一個(gè)像素位置再畫(huà)一條平行的分割線，就可以保證

在8-鄰域下割斷。下面自編的imgetline函數(shù)可通過(guò)第二個(gè)輸入?yún)?shù)n實(shí)現(xiàn)4-鄰

域或8-鄰域割斷的連續(xù)分割，n=4時(shí)為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函數(shù)，也可以用下面自編

的imdivideline函數(shù)實(shí)現(xiàn)imgetline函數(shù)的功能，兩者產(chǎn)生的分割線一般會(huì)有微小

差別。

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函數(shù)自編實(shí)現(xiàn)的4-鄰域或8-鄰域下割斷的improline函數(shù)及

程序,第二個(gè)輸入?yún)?shù)n=4時(shí)為4-鄰域割斷，其它情況(包括缺省)皆為8-鄰

域割斷，與前兩個(gè)自編函數(shù)產(chǎn)生的分割線也會(huì)有微小差別。

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

執(zhí)行結(jié)果:

圖4-13分割粘連的米粒后的圖象

4.10面積與計(jì)數(shù)

可以用bwarea函數(shù)統(tǒng)計(jì)黑白圖象中白色物體占有的總像素?cái)?shù)，用bweuler歐

拉數(shù)函數(shù)統(tǒng)計(jì)黑白圖象中孤立白色物體的個(gè)數(shù)，它有第二個(gè)可選參數(shù)，為4時(shí)

按4-鄰域割斷識(shí)別獨(dú)立的物體，為8時(shí)按8-鄰域割斷識(shí)別獨(dú)立的物體，缺省值

為8o對(duì)于圖4-12中4-鄰域割斷的圖形,按8-鄰域割斷識(shí)別得到的歐拉數(shù)為-2

(一個(gè)物體減去三個(gè)黑洞)。由此就能計(jì)算出每個(gè)白色物體占有的平均像素?cái)?shù)。

黑白圖象的歐拉數(shù)定義為：

歐拉數(shù)=白色物體的個(gè)數(shù)-白色物體內(nèi)黑色孔洞的個(gè)數(shù)

為了準(zhǔn)確檢測(cè)出白色物體的個(gè)數(shù)，事先必須填補(bǔ)白色物體內(nèi)部的黑洞，用4-

鄰域割斷后不應(yīng)再補(bǔ)洞，因?yàn)檠a(bǔ)洞會(huì)使分割線消失，而8-鄰域割斷后補(bǔ)洞不會(huì)影

響分割線。當(dāng)白色物體內(nèi)部沒(méi)有黑洞時(shí)，歐拉數(shù)就等于白色物體個(gè)數(shù)。

程序：

%檢測(cè)白色物體個(gè)數(shù)

N=bweuler(BW,4)

%檢測(cè)白色區(qū)域總的象素個(gè)數(shù)

A=bwarea(BW)

%計(jì)算每個(gè)物體的平均象素個(gè)數(shù)

a=A/N

執(zhí)行結(jié)果：

N=

70

A=

1.5447e+004

a=

220.6768

歐拉數(shù)函數(shù)的第一個(gè)輸入?yún)?shù)為填補(bǔ)過(guò)黑洞的黑白圖象矩陣，第二個(gè)參數(shù)為

鄰域選擇，可選為4和8對(duì)應(yīng)于4-鄰域和8-鄰域,缺省為8-鄰域。4-鄰域比8-

鄰域能更好區(qū)分粘連的顆粒。

[[例4-82檢測(cè)米粒圖像(文件名：rice.png)中米粒的個(gè)數(shù)及每個(gè)米粒的平均面積

(以像素為單位)。

下面程序是完全按照處理米粒的方法編寫(xiě)的程序，只是用了分區(qū)繪圖并加了

文字說(shuō)明，仍不夠完善，大家可以修改得更好。

I=imread(Tice.png');

BW=edge(I,,canny,);

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

xlab&C邊緣檢測(cè))

BW=imfill(BW,,holes,);

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

xlabelC自動(dòng)填補(bǔ)黑洞)

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

執(zhí)行結(jié)果：

N=

70

A=

1.5450e+004

a=

220.7089

執(zhí)行結(jié)果:

圖4-14求米粒圖像分析的中間過(guò)程

這是用米粒的處理程序處理的，這幅圖像有兩處需要分割，同時(shí)有一個(gè)米粒

輪廓線不封閉和兩個(gè)與邊界顆粒粘連的顆粒被去掉了，大家考慮能否修改程序，