自編第8章圖像分割

88.1(a)f(x,y)是將米粒撒在暗色絨布上拍攝的，圖(b)是該圖像的灰度直方圖。該顯然，如果選擇直方圖中兩“波峰”之間“谷底”所對應的灰度值為閾值T，那么，灰度值滿值圖像g(x,y)表示，即：gx,y

fx,yTfx,y

Interest或前景（foreground）等，其余稱為背景(background)，并約定使用灰度值1來表示物體像素，使用0來表示背景像素，但式（10.1-1）C語言編那么只需對式（8.1-1）稍作修改，就可以滿足分割結果二值圖像g(x,y)中物體和背景像素灰度值的當閾值T是一個適用于整個圖像的常數(shù)時，式（8.1-1）給出的處理稱為全局閾值分割。當閾值T在一幅圖像上處理不同像素時能改變，稱為可變閾值分割。如果像素(x,y)T取決于像素(x,y)鄰域的圖像特征(例如鄰域像素灰度均值、方差)，此時的可變閾值分割又稱為局部閾值分割然后分別在每個子圖像上再使用全局閾值分割。更進一步，如果一幅圖像中的每個像素(x,y)都要計8.2(a)顯示了一幅冰山圖像，由暗背景、冰山的明亮區(qū)域和陰影區(qū)域組成，圖(b)為其灰度直T1、T2。分割后的圖像g(x,y)由下式給出：

f(x,y)g(x,y) T1f(x,y) f(x,y)圖8.2(a)冰山圖 (c)用雙閾值將圖像分割為3個區(qū)而按式（8.1-2）T1、T2作比較，故稱雙閾值分割。以此類推，一般稱為多閾值分割，在8.1.6節(jié)中予以討論。T，就可以將該圖像分割為兩個區(qū)域。圖(b)010個灰度級的加性高斯噪聲污染后的結果，圖(e)位于兩個波峰之間的中間位置（谷底）對應的灰度級作為閾值T仍可以分割該圖像。abcdeaabde8.5a)帶有噪聲的合成圖像；(b)在[0.20.6]范圍內的灰度斜坡圖像；(c)圖(a)和圖(b)的乘積；(d)、(e)、(f)為對應的（8.1-1）TG1G2兩部分。G1T值m1和m2。T=1mm 如果TnTn1 ,停8.3（圖像灰度平均值）開始，并令T=0.52%例 clearall;closef=double(f1);count=Tnext=Delta_T=0.5;T=mean2(f);done=false;while~donecount=count+1;g=f>T;Tnext=0.5*(mean(f(g))+mean(f(~g)));done=abs(T–Tnext)<Delta_T;T=g=im2bw(f, figure figure %圖8-6(a)帶噪聲的指紋圖 （variancebetweenclasses）的概念來衡量目標和背景兩類像素之間的的可分性。如果目標和背景之,pni

i1,2,,L pi0

pi

令閾值Tk，0kL-1，按照式（8.1-1）把圖像f(xy)的像素分成C1和C2兩類，其中，C1由圖像中灰度值在[0k]范圍內取值的所有像素組成，C2則由圖像中灰度值在[k+1L-1]范圍內取值的所有像k

P2(k)P(C2)

m1kiPi|C1i ik

i i0 n j

p N （8.1-1k1Pk

j j j m2k ik1Pk

i iGGmGi

2

i

P1(k)m1(k)P2(k)m2(k) 2(k)P(k)m(k)m2P(k)m P(k)P(k)m(k)m mP(k)m(k G

BBm(k)P1(k)m1(k)i

BB) k的優(yōu)劣，Otsu使用類間方差2(k和全局方差2 B2B2G

BBGG Bkargmax2k BB)B 8.2Otsu個分量pi，灰度級i應與數(shù)據(jù)或向量的下標相對應。pni

i1,2,,L2：計算圖像灰度全局均值mGi3：計算圖像灰度全局方差 3：計算圖像灰度全局方差 i 5：for0kL-1k 計算P1(k

m(k)imP(k m(kmP(k m(k 計算2(k) G P(k)1P end11：endB最大值對應多個k值，計算這幾個k值的平均值作為k*。B2k213：計算對應最佳閾值k*的可分性度量 。2G程序運用了Matlab向量化運算技巧，以提高運算速度。function[T,SM]=%OTSUTHRESHOtsu'soptimumthresholdgivena%[T,SM]=OTSUTHRESH(H)computesanoptimumthreshold,T,intherange[01]%Otsu'smethodforagivenahistogram,%SMistheSeparabilityMeasureofClass.EqualtoEM(Effectivenessmetric)in%Normalizethehistogramtounitarea.Ifhisalready%thefollowingoperationhasnoeffect.h=h/sum(h);h= %hmustbeacolumnvectorforprocessing%Allthepossibleintensitiesrepresentedinthehistogram(256for8%(imustbeacolumnvectorforprocessingbelow.)i=(0:numel(h)-1)';%CalculatethevaluesofP1forallvaluesofk.P1=cumsum(h);%Calculatethevaluesofthemeanforallvaluesofk.m=cumsum(i.*h);%GettheglobalmeanoftheimagemG=m(end);%Calculatethebetween-class%Thedenominatoraddingtheeps,floating-pointrelativeaccuracy,toavoiddividedbyzero.sigSquared=((mG*P1-m).^2)./(P1.*(1-P1)+eps);%FindthemaximumofsigSquared.Theindexwherethemaxoccurs%theoptimumthreshold.Iftheremaybeseveralcontiguousmax%averagethemtoobtainthefinalthreshold.maxSigsq=max(sigSquared);T=mean(find(sigSquared==%Normalizedtorange[01].1issubtractedbecauseMATLABindexingstartsat%butimageintensitiesstartat0.T=(T-1)/(numel(h)-1);%SeparabilitySM=maxSigsq/(sum(((i-mG).^2).*h)+level= EM]=clearall;closeI=imread('Fig0807(a)_polymersomes.tif');[T1,EM]=graythresh(I);g1=im2bw(I,%[h,x]=[T2,SM]=otsuthresh(h);g2=im2bw(I,T2);figure,imshow(I);%顯示圖8.7(a).figure,bar(h); %顯示圖8.7(b)figure,imshow(g1);%顯示圖8.7(d).figure,imshow(g2);acac圖8.7a)聚合物囊泡的光學顯微鏡圖像；(b)直方圖；(c)使用8.1.2節(jié)的基本全局算法得到的分割結果；(d)使用Otsu圖8.8(a)是一幅簡單的合成圖像，由兩個明暗差異較大的均勻區(qū)域組成，沒有噪聲。圖(b)是它abcdefgh用Otsu方法閾值分割后的結果。clearT=graythresh(f);g=im2bw(f,%subplot(3,3,1);imshow(f); 顯示圖8.8(a).title('Originalimage');subplot(3,3,2 %顯示圖subplot(3,3,3);imshow(g); 顯示圖8.8(c).title('Segmentedresult');fn=imnoise(f,'gaussian',0,0.038);subplot(3,3,4);imshow(fn); 顯示圖8.8(d).title('Noisedimage');subplot(3,3,5);imhist(fn); 顯示圖8.8(e).Tn=graythresh(fn);gn=im2bw(fn,subplot(3,3,6);imshow(gn); 顯示圖8.8(f).title('Segmentedresult');%采用55均值濾波器對噪聲圖像進行平滑，Smooththeimage.h=fspecial('average',5);fa=imfilter(fn,h,subplot(3,3,7);imshow(fa); 顯示圖8.8(g).title('Smoothedimage');subplot(3,3,8);imhist(fa); 顯示圖8.8(h).Ta=graythresh(fa);ga=im2bw(fa,subplot(3,3,9 %顯示圖%abcde平滑后圖像的直方圖；(f)用Otsu方法對圖像進行閾值分割后的結果。兩種情形下，閾值分割均告失敗。image）用來選擇f(x,y)中參與直方圖統(tǒng)計的像素。n%對應的百分位數(shù)就是直方圖中灰度值從低到高累積概率大于n%的最小灰度值，即為要選定的閾8.7算該圖像的梯度幅度圖像，然后取其99.7%百分位數(shù)對應的梯度幅度值為閾值，對梯度幅度圖像進行分割處理，得到邊緣圖像gT(x,y)，如圖(c)所示。把圖(a)f(xy)gT(xy)相乘，結果如圖(d)所示。圖(e)是統(tǒng)計上述乘積圖像0對直方圖的貢獻。注意，該直方圖具有被一個較深的波谷分開的近乎對稱的波峰模,,abcde((方圖，用Otsu閾值分割圖像的結果。閾值為134，大約位于直方圖中兩個波峰中間位置。clearall;closeI=imread('Fig0809(a)_septagon_small_noisy_mean_0_stdv_10.tif');f=im2double(I);hr=%subplot(2,3,1);imshow(I); 顯示圖8.10(a).title('Originalimage');subplot(2,3,2);bar(hr,'hist'); 顯示圖8.10(b).%采用sobelsx=fspecial('sobel');sy=sx';gx=imfilter(f,sx,'replicate');gy=imfilter(f,sy,'replicate');grad=sqrt(gx.*gx+grad %h=h=h/sum(h); cumh=cumsum(h);row %指定直方圖99.7%百分位數(shù)對應的梯度幅度值為閾值T=row(1)/255; %此處梯度幅度已被映射到[0,255]之間，重新變換到[0,1]范圍maskImage=grad>T; %得到掩膜圖像subplot(2,3,3);imshow(maskImage); title('maskImage');fp=f.*maskImage; %屏蔽掉非邊緣像素subplot(2,3,4);imshow(fp); %顯示圖8.10(d).hp=hp(1)=0; subplot(2,3,5);bar(hp,'hist'); [Tf,SM]=otsuthresh(hp); g=im2bw(f,Tf);subplot(2,3,6imshow(g) %顯示圖%8.88.11(a)顯示了一幅酵母細胞的灰度f(xy)，希望用全局閾值分割得到圖中亮點區(qū)域。圖(b)顯示了該圖像的直方圖，圖(c)Otsu方法直接對圖像進行全局閾值分割的結果。可以看到，Otsu方法未能分割出亮點區(qū)域，盡管該方法能將一些細胞區(qū)域整體分割出來，但右側的幾個細胞區(qū)域粘連在一起沒有分開。此時，Otsu方法計算出的閾值是42，可分性度量是0.636。gT(x,y)f(x,y)與相乘，統(tǒng)計上述乘積圖像中非零像素的直方圖，如圖(e)所abcde8.11a)酵母細胞圖像；(b)圖(a)的直方圖；(c)依據(jù)圖(b)Otsu方法對圖(a)進行分割的結果；(d)取依據(jù)圖(e)中的直方圖，用Otsu方法對圖(a)進行閾值分割的結果。%8.8clearall;closeI=imread('Fig0811(a)_yeast.tif');f=im2double(I);hr=[Tr,SMr]=graythresh(f);gr=im2bw(f,Tr);subplot(2,3,1);imshow(I); 顯示圖8.11(a).title('Originalimage');subplot(2,3,2 %顯示圖subplot(2,3,3);imshow(gr); 顯示圖8.11(c).title('Segmentedresult');w=[-1-1-1;-18-1;-1-1- %Laplacianlap=abs(imfilter(f,w,'replicate'));%計算拉普拉斯絕對值圖像lap=lap/max(lap(:)); h=h=h/sum(h); cumh=cumsum(h);rowfind(cumh>=0.995);%取99.5%T=maskImagelap %fp=f.*maskImage;figure,imshow(fp);hp=imhist(fp);hp(1)=0;[Tf,SM]=otsuthresh(hp);g=im2bw(f,Tf);subplot(2,3,4);imshow(fp); 顯示圖8.11(d).subplot(2,3,5);bar(hp,'hist'); subplot(2,3,6 %顯示圖%8.1.3Otsu單閾值方法的推廣，公式中的符號含義在沒用特別指出時，同8.1.3節(jié)。Lf(xy)的灰度級數(shù)，圖像中每個像素灰度值在[0,L-1]pni

i1,2,,L我們選擇兩個閾值k1和k20<k1k2L-1，按照式（8.1-2）把圖像f(xy)中的像素劃[k1+1,k2]之間取值的像素組成，C3則由圖像中灰度值在[k2+1,L-1]之間取值的像素組成。這三類像素 G G 2Pmm2Pmm2+Pm G G

k22ik1

ik2

1 1 m1 ipi，m2 ipi，m3 1 2ik1 3ik2

（8.1- 0k1k2

加（k2k1+1，…，L-2k11，k2k1L-1的范圍內依次增加。重組或矩陣中，其行、列下標分別對應于k1和k2。g(x,y)

f(x,y) kf(x,y)k f(x,y)

最后，在2k,kk,k

GG

G圖8.12(a)冰山圖像 在MatlabR2012b及以上版本圖像處理工具箱中給出了基于Otsuthresh=multithresh(I)thresh=multithresh(I,N)[thresh,metric]=1個1N的向量，在[0,255]范圍內取值；metric為對應最佳閾值下的可分性度量。quant_I=imquantize(I,thresh);quant_I=imquantize(I,thresh,values);個區(qū)域中像素的灰度值，以式(8.1-23)給出的雙閾值分割為例，values=[c,b,a]。如果調用時忽略輸values，那么，imquantize()函數(shù)采用默認值values=[1,2,,N+1]。%8.9OtsuI=imread('Fig0812(a)_iceberg.tif');figure,imshow(I); %顯示圖8.12(a)figure,imhist(I); %顯示圖8.12(b)[thresh,metric]=multithresh(I,2);values=[0,125,255];g=imquantize(I,thresh,figure %%算法 Otsu最佳雙閾值算法Matlab實 [T1,T2,SM]=%OTSUTHRESH2Otsu'soptimumdualthresholdgivena%[T1,T2,SM]=OTSUTHRESH2(H)computestwooptimumthresholds,T1,T2,intherange[0%usingOtsu'smethodforagivenahistogram,%SMistheSeparabilityMeasureofClass.Equalto%Normalizethehistogramtounitarea.Ifhisalready%thefollowingoperationhasno%Allthepossibleintensitiesrepresentedinthehistogram(256for8bits).L=numel(h);ifL~=error('Thehistogramhmusthave256binsfor8bitsgrayimage!')h=h= %hmustbeacolumnvectorforprocessing%%(imustbeacolumnvectorforprocessingbelow.)i=(0:L-1)';%CalculatetheCumulativeProbabilityDistribution(cpd)ofhistogramh.cpd=cumsum(h);%Calculatethevaluesofthemeanforallvaluesofgraylevel.m=cumsum(i.*h);%Calculatetheimageglobalmean.mG=m(end);%Calculatethebetween-classvariance.Initializeitselementsto%Remember:0<T1<T2<L-1sigSquared=zeros(L-3,L-2);fork1=2:L-2fork2=k1+1:L-P1=P2=cpd(k2)-P1;P3=1-P1-P2;m1=m(k1)/(P1+eps); %Thedenominatoraddtheeps,toavoiddividedbyzero.m2=(m(k2)-m(k1))/(P2+eps);m3=(m(end)-%FindthemaximumofsigSquared.Theindex,(row,col),wherethemaximum%aretheoptimumthresholds,T1=row,T2=col.Theremaybeseveralcontiguousmax%Averagethemtoobtainthefinalthreshold.maxSigsq=max(sigSquared(:));[row,col]=find(sigSquared== %roundsthethresholdtothenearestgraylevel.%CalculatetheSeparabilitySM=CalculatethemaxSigsq/(sum(((i-mG).^2).*h)+f=imread('Fig0812(a)_iceberg.tif');figure,imshow(f); figure,bar(h,'hist'); [T1,T2,SM]=otsuthresh2(h);g=%%figureimshow(mat2gray(g%,,

N(s,t

fs,t

2

fs,tm

Taxy

Txyaxy y)的局部閾值Txy后，按下式判斷得到分割結果：gx,y

fx,y

式(8.1-26)和(8.1-27)用局部標準差xymxy（mG）的加權和來計算局部閾值，局限性較大。更一般意義上的局部閾值分割，常使用像素(x,y)的局部圖像特征和全局統(tǒng)計量的邏輯組合，來判斷像素(x,y)屬于背景還是目標，例如：gx,y

fx,y 8.10(下面我們采用局部閾值方法嘗試對圖(a)33Sxy計算mG，然后按式(8.1-29)局部圖像特征和全局統(tǒng)計量的邏輯組a30、b=1.5，結果如圖(d)所示。從圖中可以看到，局部閾值分割得到的目標abc圖%例 %Localthresholding.clearall;closea=b=%subplot(2,2,1);imshow(f); 顯示圖8.13(a)title('yeastimage');%%dual[thresh,metric]=multithresh(f,2);values=[0,125,255];gd=imquantize(f,thresh,subplot(2,2,2);imshow(mat2gray(gd)); 顯示圖8.13(b)title('dualthresholdingresult');%%Local%Computethelocalstandarddeviations,first.f=double(f);nhood=ones(3,3); %Bydefault,stdfiltusestheneighborhoodones(3).SIG=stdfilt(f,nhood);subplot(2,2,3 %顯示圖title('localstandard%%Computeglobalmeanoftheimage,mG.mG=mean2(f);%Obtainthesegmentedimage.g=(f>a*SIG)&(f>b*mG);subplot(2,2,4);imshow(g); 顯示圖8.13(d)title('localthresholdingresult');%

RiRi是連通域,i1,2,,j R，對所有的i和j,i jQRi i1,2,Q RjFALSE,對任何鄰接的區(qū)域Ri和Rj,i屬性是否滿足指定條件的邏輯判斷函數(shù)，比如，如果Ri中的所有像素擁有相近的灰度級，那么Q(Ri)=TURE，否則Q(Ri)=FALSE。QRiRjTRUERiRj。所謂鄰接區(qū)域，是指兩個區(qū)域在空間位置上是否鄰接，若Ri和Rj中像素的并集能形成一個連通集，那么我們說這兩個區(qū)域是鄰接的。QTRUEQ的測試結果為圖像，樹根對應于整幅圖像R，如圖8.16所示。注意，該圖中僅有子塊R4被進一步分裂。就是說，只要QRi RjTRUE，兩個鄰接區(qū)域Ri和Rj就能合并。這樣，就可以得到完整的背景( 2的整數(shù)次冪，這可通過對圖像2的整數(shù)次冪。需要明確的是，最小四分區(qū)域的尺寸與分割結果區(qū)域邊界的細節(jié)水平成反比，同時也會影響屬性邏輯判斷函數(shù)Q捕捉區(qū)域性質的有效性。8.17a)566566X射線波段圖像，要求把圍繞致密中心的0TRUE, a0<mbQ 式中，m是一個四分區(qū)域中像素的均值和標準差，ab

ab125，而標準差總是大于10。不妨令a=10，b=125。1，不滿足屬性條件的四分區(qū)域中的像素值都置為0。這樣做，能達到合并鄰接同質區(qū)域的效果。分區(qū)域尺寸（3232）則會導致分割結果的塊狀效應，區(qū)域之間的邊界輪廓將變粗糙，精細度abc8.17(a)NASA的哈勃望遠鏡在X射線波段拍攝的天鵝星座環(huán)超新星的圖像；(b)、(c)、(d)分別是將最小四分區(qū)域尺寸設為3232，1616和88像素時得到的分割結果functiong=splitmerge(f,mindim, SPLITMERGESegmentanimageusingasplit-and-merge G=SPLITMERGE(F,MINDIM,@PREDICATE)segmentsimageFby asplit-and-mergeapproachbasedonquadtree MINDIM(anonnegativeintegerpowerof2)specifiesthe dimensionofthequadtreeregions(subimages)allowed. necessary,theprogrampadstheinputimagewithzerosto nearestsquaresizethatisanintegerpowerof2. guaranteesthatthealgorithmusedinthequadtree willbeabletosplittheimagedowntoblocksofsize1-by- Theresultiscroppedbacktotheoriginalsizeofthe image.Intheoutput,G,eachconnectedregionislabeledwith different% Notethatinthefunctioncallweuse@PREDICATEforthe offun.PREDICATEisaauser-definedfunction.Itssyntax% FLAG=PREDICATE(REGION)MustreturnTRUEifthepixels REGIONsatisfythepredicatedefinedinthebodyof function;otherwise,thevalueofFLAGmustbe%% Padtheimagewithzerostothenearestsquaresizethatis integerpowerof2.Thisallowsdecompositiondowntoregions size1-by-[M,N]=size(f);f=padarray(f,[Q-M,Q-N],%PerformsplittingZ=qtdecomp(f,@split_test,mindim,%Then,performmergingbylookingateachquadregionand%allitselementsto1iftheblocksatisfiesthepredicate%infunction%First,getthesizeofthelargestblock.UsefullbecauseZ%Lmax=%Next,settheoutputimageinitiallytoallzeros.The%arrayisusedlatertoestablishconnectivity.g=zeros(size(f));MARKER=%Beginthemergingstage.form=1:Lmax[vals,r,c]=qtgetblk(f,Z,m);if~isempty(vals)%Checkthepredicateforeachoftheregionsofsizem-by-%withcoordinatesgivenbyvectorsrandc.fork=1:length(r)xlow=r(k);ylow=xhigh=xlow+m-1;yhigh=ylow+m-1;flag=fun(vals(:,:,k));ifg(xlow:xhigh,ylow:yhigh)=1;MARKER(xlow,ylow)=1;

%Finally,obtaineachconnectedregionandlabelitwith%differentintegervalueusingfunctionbwlabel.g=bwlabel(imreconstruct(MARKER,g));%Cropandexit.g=g(1:M, functionv=split_test(B,mindim,%THISFUNCTIONISPARTOFFUNCTIONSPLIT-MERGE.IT%WHETHERQUADREGIONSARESPLIT.Thefunctionreturnsin%logical1s(TRUE)fortheblocksthatshouldbesplit%logical0s(FALSE)forthosethatshould%QuadregionB,passedbyqtdecomp,isthecurrentdecomposition%theimageintokblocksofsizem-by-%kisthenumberofregionsinBatthispointintheprocedure.k=size(B,3);%Performthesplittestoneachblock.Ifthepredicate%(fun)returnsTRUE,theregionissplit,sowesetthe%elementofvtoTRUE.Else,theappropriateelementofvisset%v(1:k)=false;forI=1:kquadregion=B(:,:,ifsize(quadregion,1)<=mindimv(I)=false;flag=fun(quadregion);ifflagv(I)= functionflag=predicate(region)sd=std2(region);m=flag=(sd>10)&(m>0)&(m<%S=qtdecomp(f,@split_test,parameters)這里，f是輸入圖像，split_test是一個函數(shù)句柄，用來決定某個區(qū)域是否進行分裂，parameters是split_test函數(shù)要求的附加參數(shù)，參數(shù)之間用逗號分開。返回值S是包含四叉樹結構的稀疏矩陣，如果S(km)處的元素值非零，那么(km)是分解子塊的左上角下標，該子塊的尺寸大小等于S(km)的clearall;closefigure,imshow(f); g32=splitmerge(f,32,@predicate); figure,imshow(g32); %顯示圖8.16(b)g32=splitmerge(f,16,@predicate); g16=splitmerge(f,16,@predicate);figure %8.16g8=splitmerge(f,8,@predicate); %最小四分區(qū)域大小為88figure,imshow(g8); %顯示圖8.16(d)%abefcdghC(Mi)Mi對應的集水盆所圍繞區(qū)域像素集合，相當于該集水盆充，T[n]小于n的像素的集合，可視為水位上升到第n階段時水位之下像素的集合。即：T[n](s,t)|g(s,t) 換句話說，Cn(Mi)位于該集水盆內梯度值小于n的像素集合，顯然有：CnMiCMiT nn的增加，被水淹不減，因此，Cn-1Mi)CnMi)的子集，Cn(Mi)T[n]的子集。從以上分析可知，Cn-1(Mi)中像素所形成的一個連通分量，恰好也是T[n]中的一個連通分量。RC[n]

RC[max1] C(Mi

C[min+1]T[min+1]n步時,已經(jīng)構造出通分量qQ[n]與C[n-1]之間的關系有如下三種可能性（如圖8.21所示：qC[n1]是一個空集（8.21(a)所示qC[n1]C[n1]中的一個連通分量（8.21(b)所示qC[n1]C[n1]中多個連通分量（8.21(c)所示–nqqC[n1]Cn-1Mk)進行膨脹，要求滿足兩個條件：①膨qqC[n1]Cn-1Mi)、Cn-1Mj)兩個區(qū)域交匯聚合的那些像max+1。這樣設定可以阻止在水位n不斷升高的情況下水越過部分水壩。 8.22利用形態(tài)學膨脹構筑水壩示意圖。(a)n-1步淹沒的兩個積水盆水面區(qū)域；(b第n步時發(fā)生兩個集水盆水體transformL=watershed(A,488-連通LA相同的標記矩陣（labelmatrix），零值元與背景之間）形成分水嶺結構。二值圖像距離變換（DT，DistanceTransformofbinaryimage）是與D=bwdist(BW,distanceabcdefghi8.23(a)木銷圖像；(b)二值圖像，(c)對圖(b)求補的結果；(d)對圖(c)距離變換結果；(e)負距離圖像；(f)疊加分水線后的負距離圖像；(g疊加分水線后的二值圖像；(h將負距離圖像中對應背景區(qū)域的元素值設為負無窮(-inf)的“孤島效應”；(i)疊加新分水線后的負距離圖像%8.13clearI=imread('Fig0823(a)_wood_dowels.tif'); 顯示圖8.13(a)木銷圖像%h=I=%用OtsuIbw=figure,imshow(Ibw);%顯示圖8.23(b的二值圖像figure,imshow(~Ibw);%顯示圖8.23(c)的二值圖像Ibw1=Ibw;%%Computethedistancetransformofthecomplementofthebinaryimage.fdt=bwdist(~Ibw); %顯示圖8.23(d)距離變換結果 %顯示圖8.23(e)負距離圖像%L1=watershed(-ws1=(L1==0); fdtn1=mat2gray(-fdt);fdtn1(ws1)=0; %將分水線疊加到負距離圖像上 %顯示圖8.23(f)帶分水線的負距離圖像Ibw1(ws1)=0; %顯示圖8.23(g)%將負距離圖像圖像中屬于背景區(qū)域的元素值進強置為負無窮，-%fdtn2=-fdtn2(~Ibw)=-Inf;L2=watershed(fdtn2);ws2=(L2==fdtn3=mat2gray(-fdt);fdtn3(ws2)=0;(%直接使用分水嶺算法常會由于圖像噪聲和其他局部不規(guī)則結構的影響產生過度分割segmentationSobelg，接下來我們需要確定梯度圖像中內部標記的位置。對于解決圖8.24(a)電泳圖像的過度分割問題，所采用的內部標記應滿足：(1)被更高“海BW=imextendedmin(I,I是灰度圖像，h是一個非負的高度閾值；返回值BWI大小相等的二值圖像，其中取值為1的前景像素標記了I較深的局部極小值區(qū)域的位置。im=imexrendedmin(g,2);fim=g;fim(im)=最后兩行用中灰色斑把局部極小值區(qū)域位置疊加在原始圖像8.24a)上，如圖8.24(d)所示?？梢詄m=(Lim==0);圖8.24e)顯示了對內部標記圖像im進行距離變換的結果，圖8.24f)顯示了將分水線作為外部標記（用黑色表示，與內部標記一起疊加原始圖像中的結果。可見這些分水線恰位于由im標記的灰I2=imimposemin(I,I是灰度圖像，BWI大小相同的二值圖像，BW中的非零前景像素標值區(qū)域，來改善梯度圖像，圖8.24(g)給出了處理結果，如意圖中的灰度變化。g2=imimposemin(g,im|如圖8.24(h)所示。相對于圖8.24(c)而言分割結果得到顯著改善。abcdefgh圖8-24(a)電泳圖像；(b)梯度圖像；(c)對梯度幅值圖像進行分水嶺變換的過分割結果；(d)內部標記；(e)對內部標記距離變換結果；(f)對內部標記距離圖像分水嶺變換提取的分水線；(g)改進后的梯度幅值圖像；(h)最終分水嶺分%例8.14clearh=g=sqrt(imfilter(f,h,'replicate').^2+...imfilter(f,h','replicate').^2);wr=(L==0);f1=f1(wr)=im=imextendedmin(f,2);fim=I;fim(im)=175;fdt=bwdist(im);Lem=watershed(fdt);em=(Lem==0);fim(em)=0;g2=imimposemin(g,im|em);L2=watershed(g2);f2=I;f2(L2==0)=%256級灰度圖像而言，像素灰度值是一個標量，每個像素的灰度值都是一維坐標軸上[0,255]范圍的一個點。無論是全局閾值還是可變閾值分割，都是把一維坐標軸劃分為兩個或多RGB真彩色圖像為例，像素的屬性值對應于三維RGB顏色空間中的一個點。RGBmC。那么，這CRGB向量在顏色空間中對應的點形成了以均m為中心的“點云”，點云的形狀則與協(xié)方差矩陣C有關。TXm相似。常用的距離測度有：1 D(X,m)XmTXm xm2xm2xm2

D(X,m)XmTC1Xm DisctanceD(X,m)maxxRmR,xGmG,xBmB 選擇一種距離測度，對圖像中每一像素(xy)的屬性值執(zhí)行上述距離計算，按下式得到分割結果g(x,y)：gx,y

D(X,m)D(X,m)為2T的實心立方體。如圖8.14所示。8.14RGB顏色空間中滿足(a)歐幾里得距離(b)馬氏距離(c)T的平方比較。切比雪夫距離盡管計算簡單，但給出的相似特征點匯聚區(qū)域是一個立方gx,y

xR

TR

xG

TR

xB

x , x

x

xx x xx HSI、YCbCrI分量、Y分量為亮度分量，所以僅使用色調H與飽和度S，或色差分量Cb、Cr。(用了協(xié)方差矩陣中RGB各分量的方差信息，即協(xié)方差矩陣主對角線元素。圖(b)給出了采用歐幾里圖8.15(a)婚紗圖像 (b)歐幾里得距離分 (c)馬氏距離分 (d)切比雪夫距離分clearI=imread('Fig0815(a)_Bridewedding.jpeg');if(size(I,3)~=3)error('InputimagemustbeRGB.');%%%然后單擊鼠標右鍵，彈出菜單，選擇“CreateMask”mask=roipoly(I);%redplane=I(:,:,greenplane=I(:,:,blueplane=I(:,:,3);xR=redplane(mask);xG=greenplane(mask);xB=blueplane(mask);%ConcatenateR,G,Bcomponentarraysalongcolumn%formsampledatamatrixX=cat(2,xR,xG,xB);X=%[K,n]= %Kisthenumberofselected%Computeanunbiasedestimateofmean.misarowvectorhere.m=sum(X,1)/K;%SubtractthemeanfromeachrowofX.X=X-repmat(m,K,1);%Computeanunbiasedestimateofcovariancematrix%NotethattheproductisX'*XbecausethevectorsarerowsofX.C=(X'*X)/(K-1);%Determinethethresholdusingcovariancematrix%TeisusedforEuclidean%TmisusedforMahalanobis%TcisusedforChebyshevdisctanced=diag(C);sd=Te=3*max(sd);ifTe<10Te=Tm=ifTc<10Tc=%SegmenttheimageIusingtheparameterscalculated%ConvertItovectorformat[M,N,Q]=size(I);numpix=Ix=reshape(I,numpix,Q);f=double(Ix);fx=f-repmat(m,numpix,%ComputetheEuclideandistancebetweenallrowsofXand%DisavectorofsizeMN-by-1containingthedistance%fromallthecolorpixelstovectorm.De=sqrt(sum(fx.^2,2));%Findthedistances<=J1=find(De<=g1=zeros(M*N,1);g1(J1)=1;%Reshapeg1intoanM-by-Nimage.g1=reshape(g1,M,N);Ixm=Ix;Ixc=Ix;title('SegmentedresultusingEuclidean%%UsetheChebyshevdistancetosegmentthe%ComputetheChebyshevdistancebetweenallrowsofXandm.Dc=max(abs(fx),[],2);J3=find(Dc<=Tc);g3=zeros(M*N,1);g3(J3)=1;g3=reshape(g3,M,N);title('SegmentedresultusingChebyshev%%UsetheMahalanobisdistancetosegmentthe%ComputetheMahalanobisdistancebetweenallrowsofXandm.if(det(C)==0)error('theinversematrixofCisnotexist!');%Dm=sqrt(sum(fx*(C)^(-1).*fx,2));J2=find(Dm<=Tm);g2=zeros(M*N,1);g2(J2)=1;g2=reshape(g2,M,N);Im=reshape(Ixm,M,N,3);title('SegmentedresultusingMahalanobis%8.12Detection背景圖像減除法首先要建立描述場景圖像特征的背景模型（backgroundmodel），然后將采集的f2(xy)，…，ft(xy)fr(xy)ft(xy)進行運動目標分割，得到的結果可用二值圖像gt(x,y)表示，即：tgx,yt

ftx,yfrx,y

T2561525f(xy)為彩色圖McFarlane提出了一種基于近似中值濾波器（ApproximateMedianFilter，又稱滑動中值濾波RunningmedianFilter）ft(xy)表示像素(x,y)t的灰度值，Bt(x,y)0t

Btx,y

Bx,yBx,y

fx,yBx,y Btx,y ftx,yBtx,y 運動速度來在01之間取值。對t時刻的圖像ft(x,y)按下式將每個像素進行分類，運動目標分割結果 gx,y

ftx,yBtx,y

高斯混合模型（GMM，GaussianMixtureModels（MOG，MixtureofGaussiansf1(xy)，f2(xy)，…，ft(xy，…Xtft(xy)中任一像素(xy)的屬性值，如像素灰度值、RGB分量（rgb、HSI）等，XtRGB3Xt=[xRt,xGt,xBt]T。由于場景的變化，Xt通常Kp(Xt)wi,tXtμi,t,Σi,ti1K

K 0wi,t1,wi,t μ 1 T μXtμi,t,Σi,t

e2

i,t

i,t 3

1i,t

2I

0 22 22

i,t

0i,t

0

00 接下來，用每個像素屬性值的觀測歷史數(shù)據(jù)序列{X1，X2，…，Xt}K個和協(xié)方差矩陣Σi一定程度上揭示了第i個高斯分量屬于背景過程的可能性大小。Stauffer用各高斯分量的權重系數(shù)與標準差w/定過程的高斯分量子集。B由下式確定：BargminbwT

k 式中，TB個高斯分若T值選的較大，則背景模型往往為多峰分布，適用于背景較復雜的情形。K個高斯分量Xtk2.5Xt與該高StaufferXtXt8.4節(jié)中距離測度來XtμiXt就與該高斯分量匹配。此屬于前B個高斯分量，或不存在匹配分量，那么判定該像素為運動目標。量的參數(shù)應實時更新。標準EM算法（ExpectationMaximization，期望最大化）為批處理算法，不適合實時流數(shù)據(jù)。Stauffer采用在線K-Means聚類（onlineK-means）結合無限脈沖響應濾波器，對模型控制變量：K、α、T、V0；wi,01μi,0xR0 xB0 B=

斯分量后停止，記錄每個高斯分量的匹配結果Mk,t。Mk,t

(Xtμ)

2

gx,y wk,t1wk,t1μk,t1μk,t1 (1 (1 x G,k,t G,k,tR,k RG,k,t G,k,tG,kG,k

x B,k B,k,t B,k B,k,t B,k,t

x rateXtμk,t,Σk,t μi,tμi,ti,ti,t

wK,tμK,t R,K G,K B,K wi ,i1,,

wm值系數(shù)與標準差比值w/由大到小排序。如果Xt各分量相互獨立、但方差不相等，那么可按 BargminbwT k （framesecondabc結果；(d)采用混合高斯背景模型分割結果（運動目標用黑色顯示）abc(%Thism-fileimplementsthebackgroundsubtractionusingreference%formovingobjectsegmentation.clearall;close%Constructavideoreaderclasstoreadaavifile,firstthe'Video0801_car_parking.avi'%thenthevideoObj=VideoReader('Video0801_car_parking.avi');numFrames=videoObj.NumberOfFrames;%Readthefirstframeinthevideosequenceasthereferencebackgroundimagenewframe=read(videoObj,1);Iref=%Gettheheight,width,andnumberofcolorcomponentsoftheframe[height,width,numColor]=size(newframe);%AssignavaluetothethresholdThreh=20;fg=zeros(height,%Toavoidconsumingtoomuchmemories,readonlyaoneframeeachtime.forn=1:numFramesnewframe=read(videoObj,%Calculatetheabsolutedifferrenceimagebetweenthenew%andthereferrenceframeIdiff=abs(double(newframe)-%motionsegment,detectionmovingobjectbythreholdingIdifffg=Idiff>Threh;if(numColor== %colorfg=fg(:,:,1)|fg(:,:,2)|fg(:,:,

subplot(1,2,1),imshow(newframe);title(strcat('CurrentImage,FrameNo.',int2str(n)));subplot(1,2,2),imshow(fg);title('Segmentedresultusingreference %ApproximateMedianFilterbackground%Thism-fileimplementsthebackgroundsubtractionusingreference%formovingobjectsegmentation.clearall;close%Constructavideoreaderclasstoreadaavifile,firstthe'Video0801_car_parking.avi'%thenthevideoObj=VideoReader('Video0802_highwayII_raw.avi');numFrames=videoObj.NumberOfFrames;FPS= %GetthespeedoftheAVImovieinframespersecond%Readthefirstframeinthevideosequenceastheinitialvaluenewframe=read(videoObj,1);fmed=%Gettheheight,width,andnumberofcolorcomponentsoftheframe[height,width,numColor]=size(newframe);%AssignavaluetothethresholdThreh=20;beta=fg=zeros(height,%Toavoidconsumingtoomuchmemories,readonlyaoneframeeachtime.forn=1:numFramesnewframe=read(videoObj,%Calculatethedifferrenceimagebetweenthenew%andthereferrenceframeIref.Idiff=double(newframe)-fmed;%UpdatethemedianofeachpixelvaluepixInc=find(Idiff>0);fmed(pixInc)=fmed(pixInc)+beta;pixDec=find(Idiff<0);fmed(pixDec)=fmed(pixDec)-beta;%Motionsegment,detectionmovingobjectbythreholdingIdifffg=abs(Idiff)>Threh;if(numColor== %colorfg=fg(:,:,1)|fg(:,:,2)|fg(:,:,

subplot(1,2,1),imshow(newframe);title(strcat('CurrentImage, FrameNo.',int2str(n)));subplot(1,2,2),imshow(fg);title('SegmentedresultusingApproximateMedian %Thism-fileimplementsthemixtureofGaussiansalgorithmfor%subtractionbyusinggrayclearall%ConstructavideoreaderclasstoreadaavifilevideoObj=VideoReader('Video0802_highwayII_raw.avi');numFrames=videoObj.NumberOfFrames;%GetthespeedoftheAVImovieinframespersecond(fps)FPS=videoObj.FrameRate; framesizenewframe= %readin1stfr_gray=rgb2gray(newframe); %convertbackgroundtogreyscaleframe_size=size(fr_gray);height=frame_size(1); =frame_size(2); mog%controlK=3; %numberofgaussiancomponents(typically3-5)D=2.5; %PositivedeviationthresholdDp2= %Speedupthealpha=0.01; %learningrate(between0and1)(frompaper0.01)Thresh=0.5; %foregroundthreshold(0.25or0.75inpaper)var_init= %initialstandarddeviation(fornewcomponents)var=%Initializecomponentmeans,varianceandw=(1/K)*ones(height,width,K); %initializeweightsarrayfmean=255*ones(height,width,K); %pixelmeansfmean(:,:,1)=fvariance=var_init* %pixelu_diff=zeros(height,width,K); %differenceofeachpixelfrommeanu_diffp2=zeros(height,width,K); %savethe(u_diff*u_diff)tospeeduppho=alpha/(1/K); %initialphovariable(usedtoupdatemeanandvariance)B=3*ones(height,width); %numberofbackgroundcomponentsprocessframesforn=1:numFramesnewframe=read(videoObj, %readanew = %convertframetofg=zeros(height,%calculatedifferenceofpixelvaluesfromforu_diff(:,:,m)=abs(fr_gray-fmean(:,:,u_diffp2 %processingeachpixelinthenewframefori=1:heightforj=1:%%Gaussiancomponentmatchcheckingmatched=zeros(1,K);forc=1:if(u_diffp2(i,j,c)<=Dp2*fvariance(i,j,c)) %pixelmatchescomponentmatched(c)=1; %%%Classifythepixelif(sum(matched)==0||c>B(i,j))%noneoftheKGaussianmatchthecurrentpixel%orthematchedGaussiannotbelongtothefirst%distributions,thepixelclassifiedasforeground,andassign1fg(i,j)=1;%%Updateweights,mean,variance, %findamatchforc=1:ifw(i,j,c)=(1-alpha)*w(i,j,c)+pho %pho=alpha*(1/sqrt(2*pi*fvariance(i,j,c))) *exp(-0.5*fmean(i,j,c)=(1-pho)*fmean(i,j,c)+pho*fr_gray(i,j);fvariance(i,j,c)=(1-pho)*fvariance(i,j,c)+pho*u_