【正文】 的平行平面將生物體中需要研究的部位切成簿薄的一片，每一片就是生物體某一橫斷面的圖像，當(dāng)人們需要了解生物體信息時(shí)，再采用合成的方法，利用切片信息重建原物體的三維形態(tài)。本文為就在此背景下提出下面的問題。例如，將樣本染色后切成厚約1的切片，在顯微鏡下觀察該橫斷面的組織形態(tài)結(jié)構(gòu)。根據(jù)拍照并采樣得到的平行切片數(shù)字圖象，運(yùn)用計(jì)算機(jī)可重建組織、器官等準(zhǔn)確的三維形態(tài)。例如圓柱就是這樣一種管道，其中軸線為直線，由半徑固定的球滾動(dòng)包絡(luò)形成。圖象文件名依次為、…、格式均為，寬、高均為512個(gè)象素（）。取坐標(biāo)系的軸垂直于切片，第1張切片為平面，第100張切片為平面。問題二：確定管道的中軸線方程及畫出中軸線在平面的投影圖及三維立體圖形。二、模型假設(shè)結(jié)合本題的實(shí)際，為了確保模型求解的準(zhǔn)確性和合理性，我們排除了一些因素的干擾，提出以下幾點(diǎn)假設(shè)：所有數(shù)據(jù)均是準(zhǔn)確的，根據(jù)像素能夠近似地描繪出圖形；切片的厚度為一個(gè)像素；將血管看作一類特殊的管道，不考慮血管的彈性，即血管的半徑為一常數(shù)；對(duì)切片拍照的過程中不存在誤差，數(shù)據(jù)誤差僅與切片數(shù)字圖象的分辨率有關(guān)；中軸線上任意兩點(diǎn)處的法截面圓不相交。本文的主要工作是求解管道的半徑及中軸線方程，繪制中軸線在平面的投影圖并重建血管的三維圖像。根據(jù)圖片分析，除去圖片“彎月”的兩端之外，這一圖形的寬度是保持不變的，這一寬度就是我們所要求的球的直徑。題目中只給出了100張管道的切面圖，這些二維圖形是由無窮多個(gè)球被截的圓疊加而成的，基于：1）球的任意截面都是圓；2）經(jīng)過球心的球截面是所有截圓當(dāng)中半徑最大的圓；3）管道中軸線與每張切片有且只有一個(gè)交點(diǎn)，即為球心。首先，運(yùn)用讀取圖片，將格式文件轉(zhuǎn)化為01矩陣，通過函數(shù)確定輪廓和骨架的位置，通過循環(huán)不斷搜索計(jì)算每張圖中骨架上每一點(diǎn)到輪廓上每一點(diǎn)的最短距離，然后取最短距離中的最大值，即為最大內(nèi)切圓的半徑，最后對(duì)求出的100個(gè)最大內(nèi)切圓半徑取算數(shù)平均值減小誤差，最后的值即為管道半徑。根據(jù)擬合出曲線的光滑度確定最高階次，由于一個(gè)只對(duì)應(yīng)于一個(gè)和一個(gè),故可分別對(duì)其投影在平面上進(jìn)行多項(xiàng)式擬合，求出方程和，則中軸線的空間方程即為上兩式的聯(lián)立便得到了血管管道的中軸線參數(shù)函數(shù)。對(duì)于問題三，基于問題一、二的求解結(jié)果，根據(jù)中軸線的參數(shù)方程及滾動(dòng)球的半徑，我們運(yùn)用中的繪圖工具，結(jié)合每張切片圖的球心坐標(biāo)與球的半徑，得到血管的還原三維立體圖形。算法步驟：：導(dǎo)入數(shù)據(jù)，運(yùn)用讀取函數(shù)，轉(zhuǎn)換存儲(chǔ)方式；：運(yùn)用函數(shù)確定輪廓和骨架的位置；：尋找最大內(nèi)切圓；：求解100張圖片的最大內(nèi)切圓；：取100個(gè)半徑的算術(shù)平均值。所以在每張切片上總存在著以交點(diǎn)為圓心的圓，此時(shí)半徑最大，這最大半徑同時(shí)也是管道的半徑。（1）導(dǎo)入數(shù)據(jù)，轉(zhuǎn)換存儲(chǔ)方式本題中給了100張圖片，雖然易于直接觀察，但不利于以數(shù)學(xué)思維建模，為了減少計(jì)算量，并保證數(shù)據(jù)的準(zhǔn)確性，首先對(duì)格式的圖像進(jìn)行預(yù)處理，由于格式文件在計(jì)算機(jī)中是以二進(jìn)制數(shù)進(jìn)行存儲(chǔ)的，圖像保存在一個(gè)二維的有0或1組成的矩陣中。對(duì)于一個(gè)像素的圖像，每一個(gè)像素都有自己的一個(gè)確定的坐標(biāo)。（2）確定輪廓線【1】及骨架【2】的位置為了在減少計(jì)算量的基礎(chǔ)上，卻又反應(yīng)圖像的基本性質(zhì)，我們確定切片的輪廓線及骨架的位置，其中輪廓是指在亮度不同的區(qū)域之間有一個(gè)明顯的變化，即明度級(jí)差突然變化而形成的；而圖像的骨架即圖像中央的骨骼部分，是描述圖像幾何及拓?fù)湫再|(zhì)的重要特征之一。運(yùn)用邊緣檢測(cè)大幅度地減少了數(shù)據(jù)量，并且剔除了可以認(rèn)為不相關(guān)的信息，保留了圖像重要的結(jié)構(gòu)屬性。（3）尋找最大內(nèi)切圓 求出骨架到血管輪廓邊界上每一點(diǎn)的距離, 則其中的最小距離為即是以該點(diǎn)為圓心的內(nèi)切圓。若一張切片出現(xiàn)多個(gè)圓心時(shí)，鑒于所給切片是格式的圖像, 是對(duì)原血管的切片的連續(xù)圖形離散化而得到的近似圖形, 則實(shí)際計(jì)算時(shí)有可能出現(xiàn)在一個(gè)切片中同時(shí)找到多個(gè)最大內(nèi)切圓, 即內(nèi)切圓不惟一。在矩陣內(nèi)找到輪廓的內(nèi)切圓記下圓心坐標(biāo)與半徑并儲(chǔ)存。最后取最小距離的最大值，即是最大內(nèi)切圓的半徑：（4）100張圖片的半徑完成骨架上任一點(diǎn)對(duì)輪廓線像素點(diǎn)的遍歷求距離, 并比較記錄的所有最小距離值, 最小距離的最大值即為該切片圖像的最大內(nèi)切圓的半徑，通過循環(huán)100次可以得到這100張圖片的最大內(nèi)切圓的半徑。前兩種擬合方式顯然與現(xiàn)實(shí)相差太遠(yuǎn)，而第四種擬合又可以用第3種來逼近，所以選擇第3種方式。根據(jù)這兩個(gè)原則確定方程的最高階次，此時(shí)擬合出的中軸線較為光滑。由此畫出中軸線在平面投影擬合曲線及立體圖形，見下圖6： 圖3 中軸線在平面的投影擬合曲線 圖4中軸線在XY平面的投影擬合曲線 圖5 中軸線在平面的投影擬合曲線 圖6中軸線擬合的立體圖形分析圖像可知：1)通過三張擬合圖發(fā)現(xiàn)，得到的曲線都基本經(jīng)過散點(diǎn)圖中各三點(diǎn)的位置，特別在開始的區(qū)域，可以說是完全重合，從圖案中發(fā)現(xiàn)擬合度較高。結(jié)合圖3和圖5分析隨著的不斷變大，圖形盤旋上升的趨勢(shì)為：平—陡—平。3）盡管所得的差值沒有具體的幾何意義，但是如果對(duì)不同層的差值進(jìn)行比較，可發(fā)現(xiàn)中軸線的擬合在層數(shù)的地方誤差較小，隨著值的增大，誤差變大，這在直觀上是可以理解的；當(dāng)增大是，軸線與平面接近平行，變大，切片上容納的信息有限，產(chǎn)生較大誤差。用原來100張的平面截新的血管，生成新的100張橫斷面圖像，抽樣比較新、舊100張橫截面圖像之間的差別。所以我們可以得出結(jié)論：本模型是可行的，結(jié)果較為準(zhǔn)確。模型不能夠?qū)Π霃阶兓那闆r適應(yīng)，這是因?yàn)槲覀兊睦碚摳鶕?jù)就是等徑管道的幾何特性，另外，由于實(shí)際中圖像邊界上的點(diǎn)是連續(xù)的，只有位置，沒有大小，圖像在轉(zhuǎn)換為格式時(shí)，像素表示的邊界是離散的，產(chǎn)生了系統(tǒng)誤差。在建模的過程中，我們可以發(fā)現(xiàn)如果切片的厚度可以進(jìn)一步縮小，我們得到的擬合曲線的精度更高，反之，對(duì)于厚度過大的切片，編程的運(yùn)行所需的時(shí)間會(huì)很長，該算法不太適用。八、模型推廣如果考慮中軸線與切片的交點(diǎn)不止一個(gè)，比如有兩個(gè)交點(diǎn)，仍然可以利用我們的離散模型去求解，因?yàn)椴缓S心的切片中一定不含有最大圓。問題本身是連續(xù)的，用離散方法去模擬，不可避免地會(huì)出現(xiàn)誤差。此模型具有一定的推廣意義，尤其是在醫(yī)學(xué)觀測(cè)領(lǐng)域，醫(yī)生在不能破壞人體構(gòu)造前提下,只能通過X 光技術(shù)獲得斷層圖像信息，然而僅憑一維的斷層圖像,醫(yī)生是很難了解其內(nèi)部的復(fù)雜結(jié)構(gòu)的，因此，序列圖像的三維重現(xiàn)技術(shù)卻使斷層圖像得到立體還原，使人體內(nèi)部構(gòu)造一目了然。 三維實(shí)體的重現(xiàn)是當(dāng)今前沿科技，在醫(yī)學(xué)上有著重要的應(yīng)用?；跀?shù)字虛擬人體血管的三維重建對(duì)教學(xué)和臨床的應(yīng)用都有一定實(shí)際意義。九、參考文獻(xiàn)[1] 劉惠，圖像邊緣檢測(cè)算法分析和實(shí)現(xiàn)，湖南師范大學(xué)，2001年6月；[2] 吳德，張紅云等，圖像骨架提取的應(yīng)用，2010年4月第四期；[3] 張德豐等編，高級(jí)語言編程，機(jī)械工業(yè)出版社，2010年01月。%存儲(chǔ)以后生成的結(jié)果for k=0:99 t=strcat (int2str (k),39。)。%將bmp格式轉(zhuǎn)化為01矩陣，黑色為0，白色為1for i=1:1:512。 j0(i,j)=1J0(i,j)。sobel39。 % 提取輪廓 sobel為這個(gè)函數(shù)中截取輪廓的意思gj=bwmorph(j0,39。,inf)。%找到輪廓灰色區(qū)域[a0,b0,c0]=find(gj)。 %輪廓灰色區(qū)域個(gè)數(shù)n=length(x0)。 %建立0矩陣為求內(nèi)切圓半徑cf=zeros(m,2)。 q1=b0(i)。 q2=y0(j)。 % 通過循環(huán)求各個(gè)骨架上的各個(gè)點(diǎn)到輪廓的距離 end [zx,zxxh]=min(jl(i,:))。 cf(i,2)=zxxh。 %找到其中最大的半徑，并把半徑和圓心坐標(biāo)存儲(chǔ)%%%%%%%%%%%%%%%%%%%%%%輸出結(jié)果x=a0(zdxh)256。%與題目所給坐標(biāo)軸對(duì)應(yīng)jieguo(k+1,1)=x。%Y軸坐標(biāo)jieguo(k+1,3)=k+1。%半徑endxlswrite(39。,jieguo)%導(dǎo)入excle方便使用數(shù)據(jù)附錄二：最大內(nèi)切圓半徑及圓心坐標(biāo)最大內(nèi)切圓圓心坐標(biāo)最大內(nèi)切圓圓心坐標(biāo)切片序號(hào)最大內(nèi)切圓半徑橫坐標(biāo)縱坐標(biāo)豎坐標(biāo)切片序號(hào)最大內(nèi)切圓半徑橫坐標(biāo)縱坐標(biāo)豎坐標(biāo)4015063417046163714114966427146163724214868437246163734314868447365158744414378457468157754513788467