讓咱們經過如下簡單步驟開始咱們的配方:前端
1.經過讀取外部的體數據文件,並經過該加載數據集數據轉換成一個OpenGL紋理。也使硬件的mipmap生成。一般狀況下,從使用一個橫截面中得到的體積數據文件存儲密度影像學檢查方法,如CT或MRI掃描。每一個CT/ MRI掃描是一個二維切片。咱們在Z方向上積累簡單的2D紋理的數組得到3D紋理。密度存儲不一樣的材料的類型,例如值在0到20之間的是空氣。當咱們有一個8位無符號數據集,咱們把數據集存儲到GLubyte類型的本地數組。若是咱們有一個16位無符號的數據集咱們能夠將其存儲到GLushort類型的本地數組。3D紋理的狀況下,除了S和T的參數,咱們有一個額外的參數R控制咱們在3D紋理下的切片。ios
std::ifstream infile(volume_file.c_str(), std::ios_base::binary); if(infile.good()) { GLubyte* pData = new GLubyte[XDIM*YDIM*ZDIM]; infile.read(reinterpret_cast<char*>(pData), XDIM*YDIM*ZDIM*sizeof(GLubyte)); infile.close(); glGenTextures(1, &textureID); glBindTexture(GL_TEXTURE_3D, textureID); glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_S, GL_CLAMP); glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_T, GL_CLAMP); glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_WRAP_R, GL_CLAMP); glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR); glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_BASE_LEVEL, 0); glTexParameteri(GL_TEXTURE_3D, GL_TEXTURE_MAX_LEVEL, 4); glTexImage3D(GL_TEXTURE_3D,0,GL_RED,XDIM,YDIM,ZDIM,0,GL_RED,GL_ UNSIGNED_BYTE,pData); glGenerateMipmap(GL_TEXTURE_3D); return true; } else { return false; }
3D紋理的過濾參數相似於咱們以前看到的2D紋理參數。Mipmap是被用於細節功能的紋理的下采樣版本的集合。若是觀衆距離被應用紋理物體很是遠,它們幫助使用下采樣紋理。這有助於改善應用程序的性能。咱們必須指定最大數目的層級(GL_TEXTURE_MAX_LEVEL),便是給定紋理生成的最大mipmap數。另外,基本層(GL_TEXTURE_BASE_LEVEL)表示當對象最近時第一級所使用的mipmap。後端
glGenerateMipMap函數的經過在以前的層級過濾減小操做生成派生數組做用。因此讓咱們說咱們有三個mipmap層級,咱們的3D紋理在層級0有256*256*256的分辨率。在層級1mipmap,0級數據經過過濾減小到一半大小至128*128*128.對於第二層級mipmap,層級一被過濾並減小到64*64*64.最後,做爲第三層級,層級2將被過濾並減小到32*32*32.數組
2.設置一個頂點數組對象和頂點緩衝區對象存儲的幾何代理切片。確保緩衝區對象使用被指定爲GL_DYNAMIC_DRAW。初始glBufferData 調用分配GPU內存以得到最大切片數。vTextureSlices數組是全局定義,它存儲由頂點紋理切片操做三角產生。glBufferData用0初始化,該數據將在運行時被動態的填充。函數
const int MAX_SLICES = 512; glm::vec3 vTextureSlices[MAX_SLICES*12]; glGenVertexArrays(1, &volumeVAO); glGenBuffers(1, &volumeVBO); glBindVertexArray(volumeVAO); glBindBuffer (GL_ARRAY_BUFFER, volumeVBO); glBufferData (GL_ARRAY_BUFFER, sizeof(vTextureSlices), 0, GL_ DYNAMIC_DRAW); glEnableVertexAttribArray(0); glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE,0,0); glBindVertexArray(0);
3.經過找到一個垂直於觀察方向的具備代理切片的單位立方體交叉點實現體的切片。這是SliceVolume函數的功能。性能
由於咱們的數據在全部三個軸具備相等大小即256*256*256.咱們使用單位立方體。若是咱們有不等尺寸的數據集咱們能夠適當的縮放單位立方體。spa
//determine max and min distances glm::vec3 vecStart[12]; glm::vec3 vecDir[12]; float
float float float float lambda[12]; lambda_inc[12]; denom = 0; plane_dist = min_dist; plane_dist_inc = (max_dist-min_dist)/float(num_slices); //determine vecStart and vecDir values glm::vec3 intersection[6]; float dL[12]; for(int i=num_slices-1;i>=0;i--) { for(int e = 0; e < 12; e++) { dL[e] = lambda[e] + i*lambda_inc[e]; } if ((dL[0] >= 0.0) && (dL[0] < 1.0)) intersection[0] = vecStart[0] + dL[0]*vecDir[0]; { } //like wise for all intersection points int indices[]={0,1,2, 0,2,3, 0,3,4, 0,4,5}; for(int i=0;i<12;i++) vTextureSlices[count++]=intersection[indices[i]]; } //update buffer object glBindBuffer(GL_ARRAY_BUFFER, volumeVBO); glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(vTextureSlices), &(vTextureSlices[0].x));
4.在渲染函數,設置混合綁定體頂點數組對象,結合shader,而後調用glDrawArrays函數代理
glEnable(GL_BLEND); glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); glBindVertexArray(volumeVAO); shader.Use(); glUniformMatrix4fv(shader("MVP"), 1, GL_FALSE, glm::value_ ptr(MVP)); glDrawArrays(GL_TRIANGLES, 0, sizeof(vTextureSlices)/ sizeof(vTextureSlices[0])); shader.UnUse(); glDisable(GL_BLEND);
它是如何工做的:code
使用3D紋理切片體繪製接近由alpha混合紋理切片總體化的體繪製。第一步是加載並經過體數據生成3D紋理。裝載體數據集以後,體切片由代理切片帶出。這些體切片被定向在垂直於觀察方向的方向。此外,咱們須要找到代理多邊形和單位立方體邊界的交叉點。這由SliceVolume函數實現。須要注意的是切片只當視圖轉動時實現。orm
咱們首先得到視線方向矢量(viewDir),它是模型試圖矩陣的第三列。模型視圖矩陣的第一列存放右向量,第二列存放上向量。咱們如今詳細介紹SliceVolume函數內部如何工做。咱們經過計算在觀看方向上8個單位頂點的最小最大距離找到當前觀看方向最小最大頂點。這些距離經過每一個單位立方體頂點與視線方向向量點積獲得。
float max_dist = glm::dot(viewDir, vertexList[0]); float min_dist = max_dist; int max_index = 0; int count = 0; for(int i=1;i<8;i++) { float dist = glm::dot(viewDir, vertexList[i]); if(dist > max_dist) { max_dist = dist; max_index = i; } if(dist<min_dist) min_dist = dist; } int max_dim = FindAbsMax(viewDir); min_dist -= EPSILON; max_dist += EPSILON;
從最近頂點走到最遠頂點,從相機只有三個惟一路徑。咱們存儲每一個頂點的全部可能路徑到一個邊表,定義以下:
int edgeList[8][12]={{0,1,5,6, 4,8,11,9, 3,7,2,10 }, //v0 is front {0,4,3,11, 1,2,6,7,5,9,8,10 }, //v1 is front {1,5,0,8,2,3,7,4,6,10,9,11}, //v2 is front { 7,11,10,8, 2,6,1,9,3,0,4,5 }, // v3 is front { 8,5,9,1,11,10,7,6, 4,3,0,2 }, // v4 is front { 9,6,10,2, 8,11,4,7, 5,0,1,3 }, // v5 is front { 9,8,5,4,6,1,2,0,10,7,11,3}, // v6 is front { 10,9,6,5, 7,2,3,1,11,4,8,0 } // v7 is front
接下來,平面交叉口爲單位立方體的12個邊索引估計距離:
glm::vec3 vecStart[12]; glm::vec3 vecDir[12]; float lambda[12]; float lambda_inc[12]; float denom = 0; float plane_dist = min_dist; float plane_dist_inc = (max_dist-min_dist)/float(num_slices); for(int i=0;i<12;i++) { vecStart[i]=vertexList[edges[edgeList[max_index][i]][0]]; vecDir[i]=vertexList[edges[edgeList[max_index][i]][1]]- vecStart[i]; denom = glm::dot(vecDir[i], viewDir); if (1.0 + denom != 1.0) { lambda_inc[i] = plane_dist_inc/denom; lambda[i]=(plane_dist-glm::dot(vecStart[i],viewDir))/denom; } else { lambda[i] = -1.0; lambda_inc[i] = 0.0; } }
最後,經過在視線方向從後端到前端移動帶出內插交點和單位立方體邊。代理切片生成後,頂點緩衝對象用新的數據更新。
for(int i=num_slices-1;i>=0;i--) { for(int e = 0; e < 12; e++) { dL[e] = lambda[e] + i*lambda_inc[e]; } if ((dL[0] >= 0.0) && (dL[0] < 1.0)) { intersection[0] = vecStart[0] + dL[0]*vecDir[0]; } else if ((dL[1] >= 0.0) && (dL[1] < 1.0)) { intersection[0] = vecStart[1] + dL[1]*vecDir[1]; } else if ((dL[3] >= 0.0) && (dL[3] < 1.0)) { intersection[0] = vecStart[3] + dL[3]*vecDir[3]; } else continue; if ((dL[2] >= 0.0) && (dL[2] < 1.0)){ intersection[1] = vecStart[2] + dL[2]*vecDir[2]; } else if ((dL[0] >= 0.0) && (dL[0] < 1.0)){ intersection[1] = vecStart[0] + dL[0]*vecDir[0]; } else if ((dL[1] >= 0.0) && (dL[1] < 1.0)){ intersection[1] = vecStart[1] + dL[1]*vecDir[1]; } else {intersection[1] = vecStart[3] + dL[3]*vecDir[3]; } //similarly for others edges unitl intersection[5] int indices[]={0,1,2, 0,2,3, 0,3,4, 0,4,5}; for(int i=0;i<12;i++) vTextureSlices[count++]=intersection[indices[i]]; } glBindBuffer(GL_ARRAY_BUFFER, volumeVBO); glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(vTextureSlices), &(vTextureSlices[0].x));
在渲染函數,合適的shader被綁定。頂點shader經過對象空間頂點位置(vPosition)乘以混合模型視圖投影(MVP)矩陣,計算出夾子空間位置。它還計算用於3D紋理座標的體數據。由於咱們渲染一個單位立方體,最小的頂點位置將是(-0.5,-0.5,-0.5)而最大頂點位置將是(0.5,0.5,0.5)。因爲咱們3D紋理查詢須要從(0,0,0)到(1,1,1)的座標。咱們添加(0.5,0.5,0.5)到對象空間頂點位置來得到正確的3D紋理座標。
smooth out vec3 vUV; void main() { gl_Position = MVP*vec4(vVertex.xyz,1); vUV = vVertex + vec3(0.5); }
而後片斷shader使用3D紋理座標進行體數據採樣(其如今爲3D紋理採用一個新的採樣類型sampler3D)來顯示密度。在建立3D紋理時,咱們指定內部格式GL_RED(glTexImage3D函數的第三個參數)。所以,咱們如今能夠經過紅色通道訪問咱們的密度。要得到灰色的shader,咱們把綠色藍色和alpha通道設爲相同的值。
smooth in vec3 vUV; uniform sampler3D volume; void main(void) { vFragColor = texture(volume, vUV).rrrr; }
在之前的OpenGL版本,咱們將體密度存儲在一個特殊的內部格式GL_INTENSITY。這在OpenGL3.3核心特性已通過時了。因此如今咱們必須使用GL_RED,GL_GREEN,GL_BLUE,或GL_ALPHA內部格式。