西南交大建模培訓(xùn)MATHEMATICAMODEL 2

1、-2- 斷面可用于了解生物組織、器官等的形態(tài)。例如，將斷面可用于了解生物組織、器官等的形態(tài)。例如，將樣本染色后切成厚約樣本染色后切成厚約1 1m mm m的切片，在顯微鏡下觀察該橫斷的切片，在顯微鏡下觀察該橫斷面的組織形態(tài)結(jié)構(gòu)。如果用切片機(jī)連續(xù)不斷地將樣本切成面的組織形態(tài)結(jié)構(gòu)。如果用切片機(jī)連續(xù)不斷地將樣本切成數(shù)十、成百的平行切片，數(shù)十、成百的平行切片， 可依次逐片觀察。根據(jù)拍照并可依次逐片觀察。根據(jù)拍照并采樣得到的平行切片數(shù)字圖象，運(yùn)用計(jì)算機(jī)可重建組織、采樣得到的平行切片數(shù)字圖象，運(yùn)用計(jì)算機(jī)可重建組織、器官等準(zhǔn)確的三維形態(tài)。器官等準(zhǔn)確的三維形態(tài)。 假設(shè)某些血管可視為一類特殊的管道，該管道的表面

2、是假設(shè)某些血管可視為一類特殊的管道，該管道的表面是由球心沿著某一曲線（稱為中軸線）的球滾動(dòng)包絡(luò)而成。由球心沿著某一曲線（稱為中軸線）的球滾動(dòng)包絡(luò)而成。例如圓柱就是這樣一種管道，其中軸線為直線，由半徑固例如圓柱就是這樣一種管道，其中軸線為直線，由半徑固定的球滾動(dòng)包絡(luò)形成。定的球滾動(dòng)包絡(luò)形成。 現(xiàn)有某管道的相繼現(xiàn)有某管道的相繼100100張平行切片圖象，記錄了管道與張平行切片圖象，記錄了管道與切片的交。切片的交。-3- 為簡(jiǎn)化，假設(shè)：管道中軸線與每張切片有且只有一個(gè)交點(diǎn)；為簡(jiǎn)化，假設(shè)：管道中軸線與每張切片有且只有一個(gè)交點(diǎn)；球半徑固定；切片間距以及圖象象素的尺寸均為球半徑固定；切片間距以及圖象象素的

3、尺寸均為1 1。 取坐標(biāo)系的取坐標(biāo)系的Z Z軸垂直于切片，第軸垂直于切片，第1 1張切片為平面張切片為平面Z=0Z=0，第第100100張切片為平面張切片為平面Z=99Z=99。Z=zZ=z切片圖象中象素的坐標(biāo)依它們?cè)谇衅瑘D象中象素的坐標(biāo)依它們?cè)谖募谐霈F(xiàn)的前后次序?yàn)槲募谐霈F(xiàn)的前后次序?yàn)椋?256, -256, z）,（-256, -255, z）,（-256，255，z），），（-255,-256,z-255,-256,z）, ,（-255,-255,z-255,-255,z）, ,（-255-255，255255，z z）（ 255,-256,z 255,-256,z）, ,（ 255,

4、-255,z 255,-255,z）, ,（255255，255255，z z）。）。 試計(jì)算管道的中軸線與半徑，給出具體的算法，并繪制中軸試計(jì)算管道的中軸線與半徑，給出具體的算法，并繪制中軸線在線在XYXY、YZYZ、ZXZX平面的投影圖。平面的投影圖。-4-Z=1Z=0Z=49Z=50Z=99Z=98-5-1 1）血管的表面是由半徑固定、球心沿著某一曲血管的表面是由半徑固定、球心沿著某一曲線（稱為中軸線）的球滾動(dòng)而形成的包絡(luò)面。線（稱為中軸線）的球滾動(dòng)而形成的包絡(luò)面。2 2）中軸線上任兩點(diǎn)處的法截面圓不相交。）中軸線上任兩點(diǎn)處的法截面圓不相交。3 3）管道中軸線與每張切片平面有且只有一個(gè)交

5、）管道中軸線與每張切片平面有且只有一個(gè)交點(diǎn)。點(diǎn)。 -6- 二值圖象的象素矩陣為二值圖象的象素矩陣為0-10-1矩陣，矩陣橫矩陣，矩陣橫縱坐標(biāo)對(duì)應(yīng)原圖象的直角坐標(biāo)系位置，縱坐標(biāo)對(duì)應(yīng)原圖象的直角坐標(biāo)系位置，0 0和和1 1分分別代表該位置上的黑或白象素。別代表該位置上的黑或白象素。 0000011001110010-7- 將管道看成是小球運(yùn)動(dòng)包絡(luò)形成的，則管將管道看成是小球運(yùn)動(dòng)包絡(luò)形成的，則管道的切片便是滾動(dòng)的小球在平面道的切片便是滾動(dòng)的小球在平面z z上的截圓形上的截圓形成的包絡(luò)。成的包絡(luò)。 -30-25-20-15-10-5051015-20-1001020小 球 沿 螺 旋小 球 沿 螺 旋

6、線 滾 動(dòng) 時(shí) ，線 滾 動(dòng) 時(shí) ，在 某 平 面 上在 某 平 面 上的 截 圓 形 成的 截 圓 形 成包絡(luò)。包絡(luò)。 -8-9-解決問(wèn)題的關(guān)鍵在于發(fā)現(xiàn)解決問(wèn)題的關(guān)鍵在于發(fā)現(xiàn)定理：在一條粗細(xì)均勻血管的任何橫截面的圖象內(nèi)，其包定理：在一條粗細(xì)均勻血管的任何橫截面的圖象內(nèi)，其包含的最大內(nèi)切圓的圓心位于中軸線上，該圓的半徑等于滾含的最大內(nèi)切圓的圓心位于中軸線上，該圓的半徑等于滾動(dòng)球的半徑。動(dòng)球的半徑?；冢夯冢? 1）球的任意截面都是圓）球的任意截面都是圓2 2）經(jīng)過(guò)球心的球截面是所有）經(jīng)過(guò)球心的球截面是所有截圓當(dāng)中半徑最大的圓截圓當(dāng)中半徑最大的圓-10-求滾動(dòng)球半徑求滾動(dòng)球半徑r r的方法的方

7、法1 1）平均法）平均法2 2）抽樣法）抽樣法3 3）極大似然法）極大似然法4 4）疊加法）疊加法求出的半徑在求出的半徑在28.528.5 30.530.5之間，為什么？之間，為什么？-11-需考慮的細(xì)節(jié)：需考慮的細(xì)節(jié)：1 1）如何讀入一個(gè)數(shù)字圖象，即如何由）如何讀入一個(gè)數(shù)字圖象，即如何由bmpbmp圖象文件圖象文件（或其它格式的圖象文件）得到象素矩陣？或其它格式的圖象文件）得到象素矩陣？ MATLABMATLAB函數(shù)函數(shù)imreadimread ( ) ( )： m=imread(m=imread(0.bmp0.bmp) )或把圖象序列轉(zhuǎn)換為三維數(shù)組或把圖象序列轉(zhuǎn)換為三維數(shù)組for b=0:

8、99for b=0:99 m(:,:,b+1)=imread(int2str(b),.bmp); m(:,:,b+1)=imread(int2str(b),.bmp);endend-12-需考慮的細(xì)節(jié)：需考慮的細(xì)節(jié)：2 2）何謂邊界點(diǎn)？）何謂邊界點(diǎn)？ 四鄰域的概念四鄰域的概念 找邊界點(diǎn)坐標(biāo)的算法找邊界點(diǎn)坐標(biāo)的算法也可用也可用MATLABMATLAB函數(shù)函數(shù): :edge()edge()如：如：m=imread(m=imread(0.bmp0.bmp) )；bw=edge(m, bw=edge(m, sobelsobel) )3 3）最大內(nèi)切圓的圓心和半徑須取小數(shù)和整數(shù)有區(qū)別）最大內(nèi)切圓的圓心和

9、半徑須取小數(shù)和整數(shù)有區(qū)別嗎？嗎？4 4）取整方法：四舍五入（）取整方法：四舍五入（round)round)，向上取整向上取整( (ceil)ceil)，向下取整向下取整( (floor)floor)-13-5 5）如何由切片圖疊合作出血管的立體圖？）如何由切片圖疊合作出血管的立體圖？方法方法1 1：利用利用plot3()plot3()clearclearfor b=0:99for b=0:99 m1=imread(int2str(b),.bmp); m1=imread(int2str(b),.bmp); m(:,:,b+1)=edge(m1,sobel); m(:,:,b+1)=edge(m1

10、,sobel);endendfor k=0:99for k=0:99 for i=1:512 for i=1:512 for j=1:512 for j=1:512 if (m(i,j,k+1)=1) if (m(i,j,k+1)=1) plot3(i,j,k+1,b-.);hold on plot3(i,j,k+1,b-.);hold onend,end,end,endend,end,end,endrotate3drotate3d，hold offhold off-14-5 5）如何由切片圖疊合作出血管的立體圖？）如何由切片圖疊合作出血管的立體圖？方法方法1 1：利用利用plot3()plo

11、t3()-15-方法方法2 2：利用利用patch()patch()for b=0:99m(:,:,b+1)=imread(int2str(b),.bmp,bmp);endm1=m(1:4:512,1:4:512,:);for i=1:size(m1,1) for j=1:size(m1,2) for k=1:size(m1,3) if m1(i,j,k)=0 m1(i, j, k)=88; else m1(i,j,k)=0;end, end, end,endms=smooth3(m1);hiso=patch(isosurface(ms,5),FaceColor,1,.75,.65,EdgeC

12、olor,none);hcap=patch(isocaps(m1,5),FaceColor,interp,EdgeColor,none);colormap(default)view(45,30), axis tightgriddaspect(1,1,.4)lightangle(45,30); lighting phongisonormals(ms, hiso)set(hcap,AmbientStrength,.6)set(hiso,SpecularColorReflectance,0,SpecularExponent,50)rotate3d-16-5 5）如何由切片圖疊合作出血管的立體圖？）如

13、何由切片圖疊合作出血管的立體圖？方法方法2 2：利用利用patch()patch()-17- 中軸線的建模：求中軸線中軸線的建模：求中軸線 與各橫斷面的交點(diǎn)與各橫斷面的交點(diǎn)和曲線擬合、逼近。參賽者使用的方法主要有：和曲線擬合、逼近。參賽者使用的方法主要有：1 1）枚舉法）枚舉法2 2）平行切線法）平行切線法3 3）外推法）外推法4 4）滾球法）滾球法5 5）投影法）投影法6 6）變換法）變換法7 7）細(xì)化法）細(xì)化法-18-方法一：以黑色點(diǎn)為圓心，方法一：以黑色點(diǎn)為圓心，以以r0r0為半徑畫圓，判斷此圓為半徑畫圓，判斷此圓周 上 的 象 素 點(diǎn) 是 否 為 黑周 上 的 象 素 點(diǎn) 是 否 為

14、黑色若全為黑色，則增大半色若全為黑色，則增大半徑再判斷，直到圓周上有白徑再判斷，直到圓周上有白色象素點(diǎn)，記錄該圓的圓心色象素點(diǎn)，記錄該圓的圓心坐標(biāo)和半徑，然后，取下一坐標(biāo)和半徑，然后，取下一個(gè)黑色象素點(diǎn)重復(fù)以上步驟，個(gè)黑色象素點(diǎn)重復(fù)以上步驟，直至黑色區(qū)域中所有象素點(diǎn)直至黑色區(qū)域中所有象素點(diǎn)都搜索完為止；所有記錄點(diǎn)都搜索完為止；所有記錄點(diǎn)中對(duì)應(yīng)于最大半徑的圓心坐中對(duì)應(yīng)于最大半徑的圓心坐標(biāo)，就是該切片的最大內(nèi)切標(biāo)，就是該切片的最大內(nèi)切圓圓心坐標(biāo)，此半徑即是血圓圓心坐標(biāo)，此半徑即是血管的半徑管的半徑 -200-1000100200-1000100200020406080100 xyz-19-方法一得

15、出的中軸線投影：方法一得出的中軸線投影：0204060801001201401601800102030405060708090100y z 0102030405060708090100-200-150-100-50050100150200z x 中軸線在中軸線在ZXZX平面投影圖平面投影圖 中軸線在中軸線在YZYZ平面投影圖平面投影圖什么原因造成內(nèi)切圓圓心的誤差如此之大？什么原因造成內(nèi)切圓圓心的誤差如此之大？ -20- 模型的檢驗(yàn)?zāi)Ｐ偷臋z驗(yàn) 不檢驗(yàn)只能說(shuō)完成問(wèn)題的一半。由于所給圖象只不檢驗(yàn)只能說(shuō)完成問(wèn)題的一半。由于所給圖象只是離散數(shù)據(jù)，以象素為單位。無(wú)論以何種方式建模，是離散數(shù)據(jù)，以象素為單位

16、。無(wú)論以何種方式建模，都是近似計(jì)算，效果如何，檢驗(yàn)很必要。還可通過(guò)檢都是近似計(jì)算，效果如何，檢驗(yàn)很必要。還可通過(guò)檢驗(yàn)，發(fā)現(xiàn)模型的誤差，修正模型，提高模型的正確性。驗(yàn)，發(fā)現(xiàn)模型的誤差，修正模型，提高模型的正確性。閱卷中發(fā)現(xiàn)的檢驗(yàn)方法主要有：閱卷中發(fā)現(xiàn)的檢驗(yàn)方法主要有：1 1）逐片比較）逐片比較2 2）法平面法）法平面法3 3）滾動(dòng)法）滾動(dòng)法-21- 模型檢驗(yàn)時(shí)也發(fā)現(xiàn)內(nèi)切圓圓心誤差較大，特別是模型檢驗(yàn)時(shí)也發(fā)現(xiàn)內(nèi)切圓圓心誤差較大，特別是后面的切片。后面的切片。 Z=30Z=0-22- 模型檢驗(yàn)時(shí)也發(fā)現(xiàn)內(nèi)切圓圓心誤差較大，特別是模型檢驗(yàn)時(shí)也發(fā)現(xiàn)內(nèi)切圓圓心誤差較大，特別是后面的切片。后面的切片。 Z=5

17、0Z=40-23- 模型檢驗(yàn)時(shí)也發(fā)現(xiàn)內(nèi)切圓圓心誤差較大，特別是模型檢驗(yàn)時(shí)也發(fā)現(xiàn)內(nèi)切圓圓心誤差較大，特別是后面的切片。后面的切片。 Z=60Z=57 -24-實(shí)際圖象邊界上的點(diǎn)是連續(xù)的，在轉(zhuǎn)換成實(shí)際圖象邊界上的點(diǎn)是連續(xù)的，在轉(zhuǎn)換成bmpbmp圖象時(shí)，象素表示的圖象邊界是離散的，成圖象時(shí)，象素表示的圖象邊界是離散的，成鋸齒狀，與實(shí)際圖象有誤差（舍入誤差）。鋸齒狀，與實(shí)際圖象有誤差（舍入誤差）。2 2）同一張切片上的最大內(nèi)切圓不唯一）同一張切片上的最大內(nèi)切圓不唯一 1 1）方法一：取平均）方法一：取平均求出同一張切片上的所有最大內(nèi)切圓的圓心，然求出同一張切片上的所有最大內(nèi)切圓的圓心，然后求平均值。后

18、求平均值。-25-1 1）方法一：取平均）方法一：取平均求出同一張切片上的所有最大內(nèi)切圓的圓心，然求出同一張切片上的所有最大內(nèi)切圓的圓心，然后求平均值。后求平均值。XY平面上的投影平面上的投影YZ平面上的投影平面上的投影ZX平面上的投影平面上的投影-26-檢驗(yàn)檢驗(yàn) Z=30Z=40Z=50Z=60-27-1021021021022),(),(1),(1NyNxNyNxyxfyxgNyxeNe-28-切片Z值32435765均方誤差0.00190.00240.00480.0040 由計(jì)算出的均方誤差情況可以看出，均方由計(jì)算出的均方誤差情況可以看出，均方誤差的平均值低于誤差的平均值低于0.5% 0.5% 。 還可以用信噪比，相異象素點(diǎn)所占的比例還可以用信噪比，相異象