位相情報を考慮したMarching Cube法†CUDAでMC法で解説したMC法のCUDA実装は,頂点の重複を許しており,データの形式としては, 頂点座標の並びでポリゴンを表している. しかし,この方式ではポリゴン間の接続などは分からない. 一般的には頂点座標と頂点の接続関係の2つに分けてメッシュデータを記述する. 例えば,以下に示す図のメッシュの場合,頂点座標を表す配列が,
となり,メッシュとしては,
となる. CUDA実装†CUDAサンプルのMC法で得られたデータを解析して形式を変換しても良いが, ここでは,CUDA実装そのものを変更してみる. 手順としては,
カーネル呼び出しのコードは, void CuMCCreateMesh(float threshold, unsigned int &nvrts, unsigned int &ntris) { // MARK:CuMCCreateMesh uint lastElement, lastScanElement; int threads = 128; dim3 grid(m_uNumVoxels/threads, 1, 1); // get around maximum grid size of 65535 in each dimension if(grid.x > 65535){ grid.y = grid.x/32768; grid.x = 32768; } // // 各ボクセルの8頂点の内外を判定し,テーブルから頂点数,ポリゴン数を求める // ClassifyVoxel2<<<grid, threads>>>(m_duVoxelCubeIdx, m_duVoxelVerts, m_duVoxelTriNum, m_duVoxelOccupied, g_dfVolume, m_u3GridSize, m_uNumVoxels, m_f3VoxelH, threshold); CUT_CHECK_ERROR("classifyVoxel2 failed"); m_uNumActiveVoxels = m_uNumVoxels; #if SKIP_EMPTY_VOXELS // 空のグリッドをスキップする場合 // グリッド内のメッシュ有無配列をScan cudppScan(m_cudppScanPlan, m_duVoxelOccupiedScan, m_duVoxelOccupied, m_uNumVoxels); // Exclusive scan (最後の要素が0番目からn-2番目までの合計になっている)なので, // Scan前配列の最後(n-1番目)の要素と合計することでグリッド数を計算 cutilSafeCall(cudaMemcpy((void*)&lastElement, (void*)(m_duVoxelOccupied+m_uNumVoxels-1), sizeof(uint), cudaMemcpyDeviceToHost)); cutilSafeCall(cudaMemcpy((void*)&lastScanElement, (void*)(m_duVoxelOccupiedScan+m_uNumVoxels-1), sizeof(uint), cudaMemcpyDeviceToHost)); m_uNumActiveVoxels = lastElement+lastScanElement; if(!m_uNumActiveVoxels){ nvrts = 0; ntris = 0; return; } CompactVoxels<<<grid, threads>>>(m_duCompactedVoxelArray, m_duVoxelOccupied, m_duVoxelOccupiedScan, m_uNumVoxels); cutilCheckMsg("CompactVoxels failed"); #endif // #if SKIP_EMPTY_VOXELS // バーテックスカウント用Scan(prefix sum)作成 cudppScan(m_cudppScanPlan, m_duVoxelVertsScan, m_duVoxelVerts, m_uNumVoxels); // 三角形メッシュ数情報をScan(prefix sum) cudppScan(m_cudppScanPlan, m_duVoxelTriNumScan, m_duVoxelTriNum, m_uNumVoxels); // Exclusive scan (最後の要素が0番目からn-2番目までの合計になっている)なので, // Scan前配列の最後(n-1番目)の要素と合計することでポリゴン数を計算 cutilSafeCall(cudaMemcpy((void*)&lastElement, (void*)(m_duVoxelTriNum+m_uNumVoxels-1), sizeof(uint), cudaMemcpyDeviceToHost)); cutilSafeCall(cudaMemcpy((void*)&lastScanElement, (void*)(m_duVoxelTriNumScan+m_uNumVoxels-1), sizeof(uint), cudaMemcpyDeviceToHost)); ntris = lastElement+lastScanElement; // // エッジごとに頂点座標を計算 // uint3 dir[3]; dir[0] = make_uint3(1, 0, 0); dir[1] = make_uint3(0, 1, 0); dir[2] = make_uint3(0, 0, 1); uint cpos = 0; for(int i = 0; i < 3; ++i){ dim3 grid2(m_uNumEdges[i]/threads, 1, 1); if(grid2.x > 65535){ grid2.y = grid2.x/32768; grid2.x = 32768; } CalVertexEdge<<<grid2, threads>>>(m_dfEdgeVrts+cpos, m_duEdgeOccupied+cpos, g_dfVolume, dir[i], m_u3EdgeSize[i], m_u3GridSize, m_uNumVoxels, m_f3VoxelH, m_f3VoxelMin, threshold); cpos += m_uNumEdges[i]; } CUT_CHECK_ERROR("calVertexEdge failed"); cutilSafeCall(cudaThreadSynchronize()); // 頂点情報を詰める cudppScan(m_cudppScanPlan2, m_duEdgeOccupiedScan, m_duEdgeOccupied, m_uNumEdges[3]); cutilSafeCall(cudaMemcpy((void *) &lastElement, (void *) (m_duEdgeOccupied+m_uNumEdges[3]-1), sizeof(uint), cudaMemcpyDeviceToHost)); cutilSafeCall(cudaMemcpy((void *) &lastScanElement, (void *) (m_duEdgeOccupiedScan+m_uNumEdges[3]-1), sizeof(uint), cudaMemcpyDeviceToHost)); m_uNumTotalVerts = lastElement + lastScanElement; nvrts = m_uNumTotalVerts; if(m_uNumTotalVerts==0){ nvrts = 0; return; } dim3 grid3(m_uNumEdges[3]/threads, 1, 1); if(grid3.x > 65535){ grid3.y = grid3.x/32768; grid3.x = 32768; } // エッジ頂点を詰める CompactEdges<<<grid3, threads>>>(m_dfCompactedEdgeVrts, m_duEdgeOccupied, m_duEdgeOccupiedScan, m_dfEdgeVrts, m_uNumEdges[3]); CUT_CHECK_ERROR("compactEdges failed"); cutilSafeCall(cudaThreadSynchronize()); // // 位相情報作成 // dim3 grid4(m_uNumActiveVoxels/NTHREADS, 1, 1); if(grid4.x > 65535){ grid4.y = grid4.x/32768; grid4.x = 32768; } uint3 numEdge = make_uint3(m_uNumEdges[0], m_uNumEdges[1], m_uNumEdges[2]); GenerateTriangles3<<<grid4, NTHREADS>>>(m_du3IndexArray, m_duVoxelTriNumScan, m_duEdgeOccupiedScan, m_duVoxelCubeIdx, m_u3EdgeSize[0], m_u3EdgeSize[1], m_u3EdgeSize[2], numEdge, m_duCompactedVoxelArray, m_u3GridSize, m_uNumActiveVoxels, m_f3VoxelH, threshold, m_uNumTotalVerts, ntris); CUT_CHECK_ERROR("GenerateTriangles3 failed"); cutilSafeCall(cudaThreadSynchronize()); } ボクセル毎の頂点数,ポリゴン数のカウントのカーネルは, /*! * ボリュームデータの参照 * @param[in] data ボリュームデータ * @param[in] p グリッドインデックス * @param[in] gridSize グリッド数 * @return ボリューム値 */ __device__ float sampleVolume2(float *data, uint3 p, uint3 gridSize) { p.x = min(p.x, gridSize.x-1); p.y = min(p.y, gridSize.y-1); p.z = min(p.z, gridSize.z-1); uint i = (p.z*gridSize.x*gridSize.y)+(p.y*gridSize.x)+p.x; return data[i]; } /*! * ボクセルごとに等値頂点位置(voxelCubeIdx),等値頂点数(voxelVerts),メッシュ数(voxelTris)を計算 * @param[out] voxelCubeIdx ボクセルの等値頂点位置(8bitマスク) * @param[out] voxelVerts ボクセルの等値頂点数 * @param[out] voxelTris ボクセルのメッシュ数 * @param[in] volume ボリュームデータ * @param[in] gridSize グリッド数 * @param[in] gridSizeShift,gridSizeMask シフトとマスク * @param[in] numVoxels 総ボクセル数 * @param[in] voxelSize ボクセル幅 * @param[in] isoValue 閾値 */ __global__ void ClassifyVoxel2(uint* voxelCubeIdx, uint *voxelVerts, uint *voxelTris, uint *voxelOccupied, float *volume, uint3 ngrids, uint nvoxels, float3 voxel_h, float threshold) { uint blockId = __mul24(blockIdx.y, gridDim.x)+blockIdx.x; uint i = __mul24(blockId, blockDim.x)+threadIdx.x; if(i < nvoxels){ uint3 gridPos = calcGridIdxU(i, ngrids); // グリッド8頂点の陰関数値を参照する float field[8]; field[0] = sampleVolume2(volume, gridPos, ngrids); field[1] = sampleVolume2(volume, gridPos+make_uint3(1, 0, 0), ngrids); field[2] = sampleVolume2(volume, gridPos+make_uint3(1, 1, 0), ngrids); field[3] = sampleVolume2(volume, gridPos+make_uint3(0, 1, 0), ngrids); field[4] = sampleVolume2(volume, gridPos+make_uint3(0, 0, 1), ngrids); field[5] = sampleVolume2(volume, gridPos+make_uint3(1, 0, 1), ngrids); field[6] = sampleVolume2(volume, gridPos+make_uint3(1, 1, 1), ngrids); field[7] = sampleVolume2(volume, gridPos+make_uint3(0, 1, 1), ngrids); // グリッド内の頂点数テーブル用のインデックス作成 uint cubeindex; cubeindex = uint(field[0] < threshold); cubeindex += uint(field[1] < threshold)*2; cubeindex += uint(field[2] < threshold)*4; cubeindex += uint(field[3] < threshold)*8; cubeindex += uint(field[4] < threshold)*16; cubeindex += uint(field[5] < threshold)*32; cubeindex += uint(field[6] < threshold)*64; cubeindex += uint(field[7] < threshold)*128; uint numVerts = tex1Dfetch(numVertsTex, cubeindex); voxelCubeIdx[i] = cubeindex; // 後の計算のためにcubeindexを記録しておく voxelVerts[i] = numVerts; // グリッド内の頂点数 voxelTris[i] = numVerts/3; // グリッド内の三角形ポリゴン数 #if SKIP_EMPTY_VOXELS voxelOccupied[i] = (numVerts > 0); // グリッド内にメッシュを含むかどうか #endif } } 陰関数値の参照にはサンプルボリュームを用いている. エッジごとに頂点座標を算出するカーネルは, /*! * エッジに沿った線形補間による陰関数値0となる頂点位置の計算 * @param[in] isolavel 閾値 * @param[in] p0,p1 両端点座標 * @param[in] f0,f1 両端点の陰関数値 * @return 頂点座標 */ __device__ float3 vertexInterp(float isolevel, float3 p0, float3 p1, float f0, float f1) { float t = (isolevel-f0)/(f1-f0); return lerp(p0, p1, t); } /*! * エッジごとに0等値面頂点を探索 * @param[out] edgeVrts 頂点列 * @param[out] edgeOccupied エッジの頂点占有配列 * @param[in] volume ボリュームデータ * @param[in] dir 探索方向 * @param[in] gridSize1 エッジ数(グリッド数+1) * @param[in] gridSize グリッド数 * @param[in] gridSizeShift,gridSizeMask シフトとマスク * @param[in] numVoxels 総ボクセル数 * @param[in] voxelSize ボクセル幅 * @param[in] isoValue 閾値 */ __global__ void CalVertexEdge(float4* edgeVrts, uint *edgeOccupied, float *volume, uint3 dir, uint3 edgeSize, uint3 ngrids, uint nvoxels, float3 voxel_h, float3 voxel_min, float threshold) { uint blockId = __mul24(blockIdx.y, gridDim.x)+blockIdx.x; uint i = __mul24(blockId, blockDim.x)+threadIdx.x; uint3 gridPos = calcGridIdxU(i, edgeSize); float3 p; p.x = voxel_min.x+(gridPos.x*voxel_h.x); p.y = voxel_min.y+(gridPos.y*voxel_h.y); p.z = voxel_min.z+(gridPos.z*voxel_h.z); // calculate cell vertex positions float3 v[2]; v[0] = p; v[1] = p+make_float3(dir.x*voxel_h.x, dir.y*voxel_h.y, dir.z*voxel_h.z); // サンプルボリュームから値を取得 float field[2]; field[0] = sampleVolume2(volume, gridPos, ngrids); field[1] = sampleVolume2(volume, gridPos+dir, ngrids); uint cubeindex; cubeindex = uint(field[0] < threshold); cubeindex += uint(field[1] < threshold)*2; if(cubeindex == 1 || cubeindex == 2){ float3 vertex, normal = make_float3(0.0); vertex = vertexInterp(threshold, v[0], v[1], field[0], field[1]); edgeVrts[i] = make_float4(vertex, 1.0f); edgeOccupied[i] = 1; } else{ edgeOccupied[i] = 0; } } 頂点情報を詰めるカーネルは, /*! * ボクセルごとにScan結果に基づいてエッジ上の頂点情報を詰める * @param[out] compactedVrts 空エッジを詰めた頂点座標配列 * @param[in] occupied エッジ占有情報(頂点が存在する:1, しない:0) * @param[in] occupiedScan エッジ占有情報から作成したPrefix Sum(Scan) * @param[in] vrts 各エッジの頂点座標(頂点が存在しない場合はすべての要素が0) * @param[in] num 総エッジ数(3x(ボクセル数+1)) */ __global__ void CompactEdges(float4 *compactedVrts, uint *occupied, uint *occupiedScan, float4 *vrts, uint num) { // インデックス uint blockId = __mul24(blockIdx.y, gridDim.x)+blockIdx.x; uint i = __mul24(blockId, blockDim.x)+threadIdx.x; if(occupied[i] && (i < num)) { compactedVrts[occupiedScan[i]] = vrts[i]; } } 最終的に頂点の接続関係を計算するカーネルは, /*! * ボクセルごとにメッシュ頂点インデックスを作成 * @param[out] voxelIdx メッシュ頂点インデックス * @param[in] voxelTrisScan ボクセルごとのメッシュ数Scan * @param[in] edgeOccupiedScan エッジごとの頂点存在Scan * @param[in] voxelCubeIdx ボクセルごとの等値面頂点存在マスク * @param[in] gridSize1 エッジ数(グリッド数+1) * @param[in] gridSize グリッド数 * @param[in] gridSizeShift,gridSizeMask シフトとマスク * @param[in] numVoxels 総ボクセル数 * @param[in] voxelSize ボクセル幅 * @param[in] isoValue 閾値 */ __global__ void GenerateTriangles3(uint3 *vertIdx, uint *voxelTrisScan, uint *edgeOccupiedScan, uint *voxelCubeIdx, uint3 edgeSizeX, uint3 edgeSizeY, uint3 edgeSizeZ, uint3 edgeNum, uint *compactedVoxelArray, uint3 ngrids, uint activeVoxels, float3 voxelSize, float isoValue, uint maxVerts, uint numMesh) { uint blockId = __mul24(blockIdx.y, gridDim.x)+blockIdx.x; uint idx = __mul24(blockId, blockDim.x)+threadIdx.x; if(idx > activeVoxels-1){ idx = activeVoxels-1; } #if SKIP_EMPTY_VOXELS uint voxel = compactedVoxelArray[idx]; #else uint voxel = idx; #endif uint3 gpos = calcGridIdxU(voxel, ngrids); uint cubeindex = voxelCubeIdx[voxel]; __shared__ uint vertlist[12*NTHREADS]; vertlist[12*threadIdx.x+0] = edgeOccupiedScan[calcGridIdx3(make_uint3(gpos.x, gpos.y, gpos.z), edgeSizeX)]; vertlist[12*threadIdx.x+1] = edgeOccupiedScan[calcGridIdx3(make_uint3(gpos.x+1, gpos.y, gpos.z), edgeSizeY)+edgeNum.x]; vertlist[12*threadIdx.x+2] = edgeOccupiedScan[calcGridIdx3(make_uint3(gpos.x, gpos.y+1, gpos.z), edgeSizeX)]; vertlist[12*threadIdx.x+3] = edgeOccupiedScan[calcGridIdx3(make_uint3(gpos.x, gpos.y, gpos.z), edgeSizeY)+edgeNum.x]; vertlist[12*threadIdx.x+4] = edgeOccupiedScan[calcGridIdx3(make_uint3(gpos.x, gpos.y, gpos.z+1), edgeSizeX)]; vertlist[12*threadIdx.x+5] = edgeOccupiedScan[calcGridIdx3(make_uint3(gpos.x+1, gpos.y, gpos.z+1), edgeSizeY)+edgeNum.x]; vertlist[12*threadIdx.x+6] = edgeOccupiedScan[calcGridIdx3(make_uint3(gpos.x, gpos.y+1, gpos.z+1), edgeSizeX)]; vertlist[12*threadIdx.x+7] = edgeOccupiedScan[calcGridIdx3(make_uint3(gpos.x, gpos.y, gpos.z+1), edgeSizeY)+edgeNum.x]; vertlist[12*threadIdx.x+8] = edgeOccupiedScan[calcGridIdx3(make_uint3(gpos.x, gpos.y, gpos.z), edgeSizeZ)+edgeNum.x+edgeNum.y]; vertlist[12*threadIdx.x+9] = edgeOccupiedScan[calcGridIdx3(make_uint3(gpos.x+1, gpos.y, gpos.z), edgeSizeZ)+edgeNum.x+edgeNum.y]; vertlist[12*threadIdx.x+10] = edgeOccupiedScan[calcGridIdx3(make_uint3(gpos.x+1, gpos.y+1, gpos.z), edgeSizeZ)+edgeNum.x+edgeNum.y]; vertlist[12*threadIdx.x+11] = edgeOccupiedScan[calcGridIdx3(make_uint3(gpos.x, gpos.y+1, gpos.z), edgeSizeZ)+edgeNum.x+edgeNum.y]; __syncthreads(); // output triangle uint numTri = tex1Dfetch(numVertsTex, cubeindex)/3; for(int i = 0; i < numTri; ++i){ uint index = voxelTrisScan[voxel]+i; uint vidx[3]; uint edge[3]; edge[0] = tex1Dfetch(triTex, (cubeindex*16)+3*i); edge[1] = tex1Dfetch(triTex, (cubeindex*16)+3*i+1); edge[2] = tex1Dfetch(triTex, (cubeindex*16)+3*i+2); vidx[0] = min(vertlist[12*threadIdx.x+edge[0]], maxVerts-1); vidx[2] = min(vertlist[12*threadIdx.x+edge[1]], maxVerts-1); vidx[1] = min(vertlist[12*threadIdx.x+edge[2]], maxVerts-1); if(index < numMesh){ vertIdx[index] = make_uint3(vidx[2], vidx[1], vidx[0]); } } } |