加入實驗室後,通過張老師的介紹,有幸與某公司合共共同完成某個項目,在此項目中我主要負責的是三維 pdf 報告生成、Dicom圖像上亮度、對比度調整以及 Dicom圖像三維重建。今天主要介紹一下完成Dicom圖像三維重建的過程以及本身的心得體會。實現Dicom三維圖像重建最主要用的VTK(Visualization Toolkit,也就是可視化工具包),因爲今天的主題不是有關VTK,因此有關VTK的學習(包括VTK介紹、使用、實列),能夠參考此連接:https://blog.csdn.net/wishchin/article/details/12996693,我的建議:先把此教程中的前3個章節看完以後,在看此教程,這樣可以更好的理解程序。接下來就讓咱們進入正題。html
VTK將在可視化過程當中常常遇到的細節屏蔽起來,並封裝了一些經常使用的可視化算法,如將面繪製中經常使用的MC(MarchingCubes)算法和體繪製中經常使用的光線投射(Ray-Casting)算法封裝成類的形式提供給使用者。這樣在進行醫學體數據的可視化時就能夠直接使用VTK中已提供的相關類git
整個項目的代碼以及掛在GitHub上:https://github.com/tgpcai/Dicom_3D_Reconstruction,以爲哈不錯的能夠給樓主點一個start~github
0.三維可視化的兩種方式算法
(1)簡單點說,三維可視化的目的就是讓人們從屏幕上看到三維圖像中有什麼東西。衆所周知,二維圖像的顯示是容易的,可是三維圖像卻不同。過去因爲技術的限制,獲得了三維圖像的數據,只能把它以二維圖像的切片的形式展現給人看,這樣顯然不夠直觀。隨着計算機科學的發展,使用計算機圖形學技術爲三維物體建模並實時渲染,實現場景漫遊變成顯示三維物體的主流方法,而過去切片顯示的方式則逐漸被邊緣化。app
(2)由計算機圖形學的知識咱們能夠知道,想顯示三維圖像中的內容,能夠對這個「內容」的表面創建一個三角形網格模型。一旦獲得了這個三角網格,那麼渲染它就可以在屏幕上看到想要的內容,同時能夠調節視角進行全方位的觀察。因此第一類三維可視化方法就是基於這種思想:首先創建網格模型,以後渲染網格。這種方式被稱爲面繪製。函數
(3)還有一種叫作體繪製的方式,是直接將三維圖像的體素點經過必定的透明度疊加計算後直接對屏幕上的像素點着色。這種方式的特色是能更加清楚的表現體數據內部細節,可是這種算法通常對計算機的壓力也會比較大。工具
1.基於面繪製的MC算法佈局
(0)首先基於MC的一系列算法須要明確一個「體元(Cell)」的概念。體元是在三維圖像中由相鄰的八個體素點組成的正方體方格,MarchingCubes算法的Cube的語義也能夠指這個體元。注意區別體元和體素,體元是8個體素構成的方格,而每一個體素(除了邊界上的以外)都爲8個體元所共享。學習
(1)面繪製:面繪製是採用分割技術對一系列的二維圖像進行輪廓識別、提取等操做,最終還原出被檢測物體的三維模型,並以表面的方式顯示出來。開發工具
(2)面繪製實現三維重建。使用的是經典的 Marching Cubes 算法,也叫移動立方體法。
(3)採用面繪製,VTK中的數據流以下:source->filter(MC算法或者vtkContourFilter)->mapper->actor->render->renderwindow->interactor。
(4)MC算法簡介:
(5)VTK提供了兩種提取等值面的類:vtkContourFilter濾波器和封裝了MC(Marching Cubes)算法類vtkMarchingCubes。提取等值面以後的數據處理:經過vtkPolyDataNormals在等值面上產生法向量;經過vtkStripper在等值面上產生紋理或三角面片。
(6)利用MC算法提取等值面的代碼實現:
1 import vtk 2 # source->filter(MC算法)->mapper->actor->render->renderwindow->interactor 3 4 # 讀取Dicom數據,對應source 5 v16 = vtk.vtkDICOMImageReader() 6 # v16.SetDirectoryName('D:/dicom_image/V') 7 v16.SetDirectoryName('D:/dicom_image/vtkDicomRender-master/sample') 8 9 # 利用封裝好的MC算法抽取等值面,對應filter 10 marchingCubes = vtk.vtkMarchingCubes() 11 marchingCubes.SetInputConnection(v16.GetOutputPort()) 12 marchingCubes.SetValue(0, 100) 13 14 # 剔除舊的或廢除的數據單元,提升繪製速度,對應filter 15 Stripper = vtk.vtkStripper() 16 Stripper.SetInputConnection(marchingCubes.GetOutputPort()) 17 18 # 創建映射,對應mapper 19 mapper = vtk.vtkPolyDataMapper() 20 # mapper.SetInputConnection(marchingCubes.GetOutputPort()) 21 mapper.SetInputConnection(Stripper.GetOutputPort()) 22 23 # 創建角色以及屬性的設置,對應actor 24 actor = vtk.vtkActor() 25 actor.SetMapper(mapper) 26 # 角色的顏色設置 27 actor.GetProperty().SetDiffuseColor(1, .94, .25) 28 # 設置高光照明係數 29 actor.GetProperty().SetSpecular(.1) 30 # 設置高光能量 31 actor.GetProperty().SetSpecularPower(100) 32 33 # 定義舞臺,也就是渲染器,對應render 34 renderer = vtk.vtkRenderer() 35 36 # 定義舞臺上的相機,對應render 37 aCamera = vtk.vtkCamera() 38 aCamera.SetViewUp(0, 0, -1) 39 aCamera.SetPosition(0, 1, 0) 40 aCamera.SetFocalPoint(0, 0, 0) 41 aCamera.ComputeViewPlaneNormal() 42 43 # 定義整個劇院(應用窗口),對應renderwindow 44 rewin = vtk.vtkRenderWindow() 45 46 # 定義與actor之間的交互,對應interactor 47 interactor = vtk.vtkRenderWindowInteractor() 48 49 # 將相機添加到舞臺renderer 50 renderer.SetActiveCamera(aCamera) 51 aCamera.Dolly(1.5) 52 53 # 設置交互方式 54 style = vtk.vtkInteractorStyleTrackballCamera() 55 interactor.SetInteractorStyle(style) 56 57 # 將舞臺添加到劇院中 58 rewin.AddRenderer(renderer) 59 interactor.SetRenderWindow(rewin) 60 61 # 將角色添加到舞臺中 62 renderer.AddActor(actor) 63 64 # 將相機的焦點移動至中央,The camera will reposition itself to view the center point of the actors, 65 # and move along its initial view plane normal 66 renderer.ResetCamera() 67 68 interactor.Initialize() 69 interactor.Start()
結果以下:
(7)利用vtkContourFilter濾波器提取等值面的代碼實現:
1 # 抽取輪廓(等值面)的操做對象是標量數據。 2 # 其思想是:將數據集中標量值等於某一指定恆量值的部分提取出來。對於3D的數據集而言,產生的是一個等值面;對於2D的數據集而言,產生的是一個等值線。 3 # 其典型的應用有氣象圖中的等溫線、地形圖中的等高線。對於醫學數據而言,不一樣的標量值表明的是人體的不一樣部分,於是能夠分別提取出人的皮膚或骨頭。 4 # 抽取輪廓的功能是由一個過濾器實現的,如vtkContourFilter、vtkMarchingCubes。vtkContourFilter能夠接受任意數據集類型做爲輸入,於是具備 通常性。 5 # 使用vtkContourFilter 時,除了須要設置輸入數據集外,還須要指定一個或多個用於抽取的標量值。可用以下兩種方法進行設置。 6 # 7 # 使用方法SetValue()逐個設置抽取值。該方法有個兩個參數:第一個參數是抽取值的索引號,表示第幾個 抽取值。索引號從0開始計數;第二個參數就是指定的抽取值。 8 # 使用方法GenerateValues()自動產生一系列抽取值。該方法有三個參數:第一個參數是抽取值的個數,後面兩個參數是抽取值的取值範圍。 9 10 # coding=utf-8 11 import vtk 12 13 # source—filter——mapper——actor——render——renderwindow——interactor 14 aRenderer = vtk.vtkRenderer() # 渲染器 15 renWin = vtk.vtkRenderWindow() # 渲染窗口,建立窗口 16 renWin.AddRenderer(aRenderer) # 渲染窗口 17 # renWin.Render() 18 iren = vtk.vtkRenderWindowInteractor() # 窗口交互 19 iren.SetRenderWindow(renWin) 20 21 # The following reader is used to read a series of 2D slices(images) 22 # that compose the volume.Theslicedimensions are set, and the 23 # pixel spacing.The data Endianness must also be specified.The reader 24 # uses the FilePrefix in combination with the slice number to construct 25 # filenames using the format FilePrefix. % d.(In this case the FilePrefix 26 # is the root name of the file. 27 28 v16 = vtk.vtkDICOMImageReader() 29 # v16.SetDirectoryName('D:/dicom_image/V') 30 v16.SetDirectoryName('D:/dicom_image/vtkDicomRender-master/sample') 31 32 33 34 # An isosurface, or contour value of 500 is known to correspond to the 35 # skin of the patient.Once generated, a vtkPolyDataNormals filter is 36 # used to create normals for smooth surface shading during rendering. 37 skinExtractor = vtk.vtkContourFilter() 38 skinExtractor.SetInputConnection(v16.GetOutputPort()) 39 skinExtractor.SetValue(0, -10) 40 # skinExtractor.GenerateValues(2, 100, 110) 41 skinNormals = vtk.vtkPolyDataNormals() 42 skinNormals.SetInputConnection(skinExtractor.GetOutputPort()) 43 skinNormals.SetFeatureAngle(60.0) 44 skinMapper = vtk.vtkPolyDataMapper() # 映射器 45 skinMapper.SetInputConnection(skinNormals.GetOutputPort()) 46 skinMapper.ScalarVisibilityOff() 47 48 skin = vtk.vtkActor() 49 # 設置顏色RGB顏色系統就是由三個顏色份量:紅色(R)、綠色(G)和藍色(B)的組合表示, 50 # 在VTK裏這三個份量的取值都是從0到1,(0, 0, 0)表示黑色,(1, 1, 1)表示白色。 51 # vtkProperty::SetColor(r,g, b)採用的就是RGB顏色系統設置顏色屬性值。 52 #skin.GetProperty().SetColor(0, 0, 1) 53 skin.SetMapper(skinMapper) 54 55 skin.GetProperty().SetDiffuseColor(1, .49, .25) 56 57 skin.GetProperty().SetSpecular(.5) 58 59 skin.GetProperty().SetSpecularPower(20) 60 61 # skin.GetProperty().SetRepresentationToSurface() 62 # 構建圖形的方框 63 outlineData = vtk.vtkOutlineFilter() 64 outlineData.SetInputConnection(v16.GetOutputPort()) 65 mapOutline = vtk.vtkPolyDataMapper() 66 mapOutline.SetInputConnection(outlineData.GetOutputPort()) 67 outline = vtk.vtkActor() 68 outline.SetMapper(mapOutline) 69 outline.GetProperty().SetColor(0, 0, 0) 70 71 # 構建舞臺的相機 72 aCamera = vtk.vtkCamera() 73 aCamera.SetViewUp(0, 0, -1) 74 aCamera.SetPosition(0, 1, 0) 75 aCamera.SetFocalPoint(0, 0, 0) 76 aCamera.ComputeViewPlaneNormal() 77 78 # Actors are added to the renderer.An initial camera view is created. 79 # The Dolly() method moves the camera towards the Focal Point, 80 # thereby enlarging the image. 81 aRenderer.AddActor(outline) 82 aRenderer.AddActor(skin) 83 aRenderer.SetActiveCamera(aCamera) 84 # 將相機的焦點移動至中央,The camera will reposition itself to view the center point of the actors, 85 # and move along its initial view plane normal 86 aRenderer.ResetCamera() 87 # aCamera.Dolly(1.5) 88 # aCamera.Roll(180) 89 # aCamera.Yaw(60) 90 91 aRenderer.SetBackground(250, 250, 250) 92 # renWin.SetSize(640, 480) 93 # 該方法是從vtkRenderWindow的父類vtkWindow繼承過來的,用於設置窗口的大小,以像素爲單位。 94 renWin.SetSize(500, 500) 95 aRenderer.ResetCameraClippingRange() 96 97 style = vtk.vtkInteractorStyleTrackballCamera() 98 iren.SetInteractorStyle(style) 99 100 iren.Initialize() 101 iren.Start()
結果以下:
(8)與可視化窗口的交互方式:可使用鼠標與三維圖形交互,好比用鼠標滾輪能夠對三維圖形放大、縮小;按下鼠標左鍵不放,而後移動鼠標,能夠轉動三維圖形;按下鼠標左鍵,同時按下Shift鍵,移動鼠標,能夠移動整個三維圖形,等等。其餘的功能你也能夠試着摸索一下,好比按下Ctrl鍵時再按鼠標左鍵;鼠標停留在柱體上,而後按下P鍵;按一下字母E將關閉窗口。
(9)整個過程的簡要總結:
(a)讀取數據以及數據處理:首先,讀取切片數據,並將其轉換爲咱們的開發工具VTK所支持的一種數據表達形式(vtkImageData)。咱們給CT數據創建的是比較抽象的等值面模型,最後將物理組件與抽象的模型結合在一塊兒來創建對CT 數據的可視化,以幫助用戶正確理解數據。利用VTK中的vtkDICOMImageReader 咱們能夠很方便的讀取切片數據。
(b)提取等值面:接着咱們就能夠用算法對所讀取的數據進行處理了。好比採用的經典MC的面繪製方法,首先利用vtkMarchingCubes 類來提取出某一CT值的等值面(利用vtksetValue()來設置須要提取的值),但這時的等值面其實仍只是一些三角面片,還必須由vtkStripper類將其拼接起來造成連續的等值面。這樣就把讀取的原始數據通過處理轉換爲應用數據,也即由原始的點陣數據轉換爲多邊形數據而後由vtkPolyDataMapper將其映射爲幾何數據,並將其屬性賦給窗口中表明它的演員,將結果顯示出來。在實際應用中Visualization Toolkit支持多表面重建。咱們能夠設置多個參數值,提取出多個等值面並同時顯示出來,如何設置多個參數值呢?能夠經過VTK自帶的GenerateValues()函數。常見的好比人體皮膚所對應的value值爲500,人體骨骼所對應的value值爲1150。
(c)顯示結果:經過前面這些工做,咱們基本上已經完成了對數據的讀取處理映射等步驟,下面咱們就要對數據進行顯示了。一般這些步驟也叫作渲染引擎。能夠經過調整value值和actor的相應屬性達到重建三維圖形的不一樣效果。
最主要就是設置相機,設置actor的相關屬性(顏色,亮度,透明度等等)。
2.基於體繪製的 Ray-casting算法
(1)體繪製:體繪製是將三維空間的離散數據直接轉換爲最後的立體,圖像而沒必要生成中間幾何圖元(面繪製須要), 其中心思想是爲每個體素指定一個不透明度,並考慮每個體素對光線的透射、發射和反射做用。
(2)體繪製達到的效果:體繪製的目標是在一副圖片上展現空間體細節。舉例而言,你面前有一間房子,房子中有傢俱、家電,站在房子外面只能看到外部形狀(相似於面繪製的效果),沒法觀察到房子的佈局或者房子中的物體;假設房子和房子中的物體都是半透明的,這樣你就能夠同時查看到全部的細節。這就是體繪製所要達到的效果。
(3)體繪製經常使用的算法:光線投射算法( Ray-casting )、錯切 - 變形算法( Shear-warp )、頻域體繪製算法( Frequency Domain )和拋雪球算法( Splatting )。其中又以光線投射算法最爲重要和通用。
(4)光線投射算法( Ray-casting )原理:從圖像平面的每一個像素都沿着視線方向發出一條射線,此射線穿過體數據集,按必定步長進行採樣,由內插計算每一個採樣點的顏色值和不透明度,而後由前向後或由後向前逐點計算累計的顏色值和不透明度值,直至光線徹底被吸取或穿過物體。該方法能很好地反映物質邊界的變化,使用Phong模型,引入鏡面反射、漫反射和環境反射能獲得很好的光照效果,在醫學上可將各組織器官的性質屬性、形狀特徵及相互之間的層次關係表現出來,從而豐富了圖像的信息。(借鑑百度百科)
(5)體繪製的原理和麪繪製徹底不相同。面繪製須要生成中間圖元,而體繪製則是直接在原圖上進行繪製,內容需求較面繪製小。每切換一個視角須要從新對全部的像素點進行顏色和透明度計算,須要時間比面繪製長。
(6)代碼實現基於體繪製的 Ray-casting算法
1 # This example reads a volume dataset and displays it via volume rendering(體繪製). 2 3 import vtk 4 from vtk.util.misc import vtkGetDataRoot 5 6 # Create the renderer, the render window, and the interactor. The renderer 7 # draws into the render window, the interactor enables mouse- and 8 # keyboard-based interaction with the scene. 9 ren = vtk.vtkRenderer() 10 renWin = vtk.vtkRenderWindow() 11 renWin.AddRenderer(ren) 12 iren = vtk.vtkRenderWindowInteractor() 13 iren.SetRenderWindow(renWin) 14 15 # The following reader is used to read a series of 2D slices (images) 16 # that compose the volume. The slice dimensions are set, and the 17 # pixel spacing. The data Endianness must also be specified. The reader 18 # usese the FilePrefix in combination with the slice number to construct 19 # filenames using the format FilePrefix.%d. (In this case the FilePrefix 20 # is the root name of the file: quarter.) 21 22 # v16 = vtk.vtkVolume16Reader() 23 # v16.SetDataDimensions(64, 64) 24 # v16.SetImageRange(1, 93) 25 # v16.SetDataByteOrderToLittleEndian() 26 # v16.SetFilePrefix("D:/dicom_image/headsq/quarter") 27 # v16.SetDataSpacing(3.2, 3.2, 1.5) 28 v16 = vtk.vtkDICOMImageReader() 29 # v16.SetDirectoryName('D:/dicom_image/vtkDicomRender-master/sample') 30 v16.SetDirectoryName('D:/dicom_image/V') 31 32 # The volume will be displayed by ray-cast alpha compositing. 33 # A ray-cast mapper is needed to do the ray-casting, and a 34 # compositing function is needed to do the compositing along the ray. 35 volumeMapper = vtk.vtkGPUVolumeRayCastMapper() 36 volumeMapper.SetInputConnection(v16.GetOutputPort()) 37 volumeMapper.SetBlendModeToComposite() 38 39 # The color transfer function maps voxel intensities to colors. 40 # It is modality-specific, and often anatomy-specific as well. 41 # The goal is to one color for flesh (between 500 and 1000) 42 # and another color for bone (1150 and over). 43 volumeColor = vtk.vtkColorTransferFunction() 44 volumeColor.AddRGBPoint(0, 0.0, 0.0, 0.0) 45 volumeColor.AddRGBPoint(500, 1.0, 0.5, 0.3) 46 volumeColor.AddRGBPoint(1000, 1.0, 0.5, 0.3) 47 volumeColor.AddRGBPoint(1150, 1.0, 1.0, 0.9) 48 49 # The opacity transfer function is used to control the opacity 50 # of different tissue types. 51 volumeScalarOpacity = vtk.vtkPiecewiseFunction() 52 volumeScalarOpacity.AddPoint(0, 0.00) 53 volumeScalarOpacity.AddPoint(500, 0.15) 54 volumeScalarOpacity.AddPoint(1000, 0.15) 55 volumeScalarOpacity.AddPoint(1150, 0.85) 56 57 # The gradient opacity function is used to decrease the opacity 58 # in the "flat" regions of the volume while maintaining the opacity 59 # at the boundaries between tissue types. The gradient is measured 60 # as the amount by which the intensity changes over unit distance. 61 # For most medical data, the unit distance is 1mm. 62 volumeGradientOpacity = vtk.vtkPiecewiseFunction() 63 volumeGradientOpacity.AddPoint(0, 0.0) 64 volumeGradientOpacity.AddPoint(90, 0.5) 65 volumeGradientOpacity.AddPoint(100, 1.0) 66 67 # The VolumeProperty attaches the color and opacity functions to the 68 # volume, and sets other volume properties. The interpolation should 69 # be set to linear to do a high-quality rendering. The ShadeOn option 70 # turns on directional lighting, which will usually enhance the 71 # appearance of the volume and make it look more "3D". However, 72 # the quality of the shading depends on how accurately the gradient 73 # of the volume can be calculated, and for noisy data the gradient 74 # estimation will be very poor. The impact of the shading can be 75 # decreased by increasing the Ambient coefficient while decreasing 76 # the Diffuse and Specular coefficient. To increase the impact 77 # of shading, decrease the Ambient and increase the Diffuse and Specular. 78 volumeProperty = vtk.vtkVolumeProperty() 79 volumeProperty.SetColor(volumeColor) 80 volumeProperty.SetScalarOpacity(volumeScalarOpacity) 81 # volumeProperty.SetGradientOpacity(volumeGradientOpacity) 82 volumeProperty.SetInterpolationTypeToLinear() 83 volumeProperty.ShadeOn() 84 volumeProperty.SetAmbient(0.9) 85 volumeProperty.SetDiffuse(0.9) 86 volumeProperty.SetSpecular(0.9) 87 88 # The vtkVolume is a vtkProp3D (like a vtkActor) and controls the position 89 # and orientation of the volume in world coordinates. 90 volume = vtk.vtkVolume() 91 volume.SetMapper(volumeMapper) 92 volume.SetProperty(volumeProperty) 93 94 # Finally, add the volume to the renderer 95 ren.AddViewProp(volume) 96 97 # Set up an initial view of the volume. The focal point will be the 98 # center of the volume, and the camera position will be 400mm to the 99 # patient's left (which is our right). 100 camera = ren.GetActiveCamera() 101 c = volume.GetCenter() 102 camera.SetFocalPoint(c[0], c[1], c[2]) 103 camera.SetPosition(c[0] + 400, c[1], c[2]) 104 camera.SetViewUp(0, 0, -1) 105 106 # Increase the size of the render window 107 renWin.SetSize(640, 480) 108 109 # Interact with the data. 110 iren.Initialize() 111 renWin.Render() 112 iren.Start()
結果以下:
(7)體繪製的整個過程包括VTK數據量都與面繪製相似。一樣能夠經過調整actor的相應屬性達到重建三維圖形的不一樣效果,好比經過設置不透明度值來顯示體數據內部的不一樣成分和細節,例如顯示人體CT圖像的不一樣器官和組織。
3.心得體會與總結
(1)其實不管是面繪製仍是體繪製都須要必定的VTK知識,因此先了解VTK的一些基礎知識才能幫助你更好的掌握這些方法。
(2)有關VTK整個數據流的過程能夠用一下的例子進行類比,方便理解(雖然這個類比不是很是形象):
(3)體繪製與面繪製管線上的區別
參考文獻:https://vtk.org/doc/nightly/html/annotated.html
http://blog.sina.com.cn/s/blog_5ff6097b0100zz2y.html
以上就是本次學習的內容,歡迎交流與討論