I=imread（‘D:\MATLAB\R2012a\bin\bmpCT庫\001.bmp‘）;
I=rgb2gray（I）;
subplot（221）;imshow（I）;
xlabel（‘a）原始圖像‘）;
%使用graythresh函數計算閾值
level=graythresh（I）;
%大津法計算全局圖像I的閾值
%BW=im2bw（Ilevel）;
%subplot（222）;imshow（BW）;
%xlabel（‘b）graythresh函數計算閾值‘）;
%disp（strcat（‘graythresh函數計算灰度閾值：‘num2str（uint8（level*255））））;
%下面matlab程序實現簡化計算閾值
IMAX=max（max（I））;
IMIN=min（min（I））;
T=double（IMIN:IMAX）;
ISIZE=size（I）;%圖像大小
muxSize=ISIZE（1）*ISIZE（2）;
for?i=1:length（T）
????%從最小灰度值到最大值分別計算方差
????TK=T（1i）;
????ifground=0;
????ibground=0;
????%定義前景和背景灰度總和
????FgroundS=0;
????BgroundS=0;
????for?j=1:ISIZE（1）
????????for?k=1:ISIZE（2）
????????????tmp=I（jk）;
????????????if（tmp>=TK）
????????????????ifground=ifground+1;
????????????????FgroundS=FgroundS+double（tmp）;%前景灰度值
????????????else
????????????????%背景像素點的計算
????????????????ibground=ibground+1;
????????????????BgroundS=BgroundS+double（tmp）;
????????????end
????????end
????end
????%計算前景和背景的比例、平均灰度值
????%這里存在一個0分母的情況，導致警告，但不影響結果
????w0=ifground/muxSize;
????w1=ibground/muxSize;
????u0=FgroundS/ifground;
????u1=BgroundS/ibground;
????T（2i）=w0*w1*（u0-u1）*（u0-u1）;%計算方差
end
%遍歷后尋找I第二行的最大值
oMax=max（T（2:））;
%第二行方差的最大值，忽略NaN
idx=find（T（2:）>=oMax）;
%方差最大值所對應的列號
T=uint8（T（1idx））;
%從第一行取出灰度值作為賦值
disp（strcat（‘簡化大津法計算灰度閾值：‘num2str（T）））;
BW=im2bw（Idouble（T）/255）;%閾值分割
subplot（223）;imshow（BW）;
xlabel（‘b）簡化大津法計算閾值‘）;
f=BW;
?F=imfill（f8‘holes‘）;%填充肺實質空洞
?E=F-f;
?E=imfill（E8‘holes‘）;%填充肺實質空洞

f=bwareaopen（E500）;%去除面積小于T部分（氣管），取值可變
[gt]=edge（f‘canny‘）;
subplot（222）;imshow（g）;
xlabel（‘c）canny‘）;
B=bwboundaries（g8‘holes‘）;
F=imfill（g8‘holes‘）;
subplot（236）;imshow（F）;

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件??????786486??2017-07-11?16:05??otsu\001.bmp
?????文件????????1930??2017-09-12?10:26??otsu\Otsu.m
?????目錄???????????0??2017-09-12?10:26??otsu\