上一篇

以xml文本的方式读取vtu文件

发布于 2023-11-23 346 次阅读

Please refresh the page if equations are not rendered correctly.
---------------------------------------------------------------

以ASCII格式存储的vtu网格文件可以采用xml文本的方式读写。然而效率比较低。当年暴力的写出来了这些代码，现在回头看不过是并不高明的造轮子。然而，这些代码对于帮助理解vtu网格文件的格式和xml文件的操作还是很有帮助的，毕竟很多文件格式本质上就是xml的方言而已。特此记录。

1. 指定单元数和尺寸创建体素网格

import numpy as np
import vtk
import xml.dom.minidom as minidom
from xml.dom import minidom

import timeit

'''
创建足够大的NumArray然后对应赋值，而不是使用np.insert()可以大大提高数据处理速度

'''
##################################
###  generate points & cells   ###
##################################

starttime1 = timeit.default_timer()

dx, dy, dz = 1, 1, 1

# number of cells in different directions
nrow = 100
ncolumn = 150
ndepth = 60

lx = dx*nrow
ly = dy*ncolumn
lz = dz * ndepth

ncells = nrow * ncolumn * ndepth
npoints = (nrow + 1) * (ncolumn + 1) * (ndepth + 1)

# Coordinates
X = np.arange(0, lx + 0.1*dx, dx, dtype='float64')
Y = np.arange(0, ly + 0.1*dy, dy, dtype='float64')
Z = np.arange(0, lz + 0.1*dz, dz, dtype='float64')

x = np.zeros((nrow + 1, ncolumn + 1, ndepth + 1))
y = np.zeros((nrow + 1, ncolumn + 1, ndepth + 1))
z = np.zeros((nrow + 1, ncolumn + 1, ndepth + 1))

for k in range(ndepth + 1):
    for j in range(ncolumn + 1):
        for i in range(nrow + 1): 
            x[i,j,k] = X[i]
            y[i,j,k] = Y[j]
            z[i,j,k] = Z[k]
# points
x = x.reshape([1,-1])[0]
y = y.reshape([1,-1])[0]
z = z.reshape([1,-1])[0]

points = np.zeros([x.size,4])
points[:,0] = x # points_organized[:,:,:,0]
points[:,1]=y # points_organized[:,:,:,1]
points[:,2]=z # points_organized[:,:,:,2]
points[:,3]=np.arange(x.size) # index of each point, points_organized[:,:,:,3]

##print('-------Points list--------')
##print(points)

# connectivity: Speed up cell allocation
points_organized = points.reshape([nrow+1,ncolumn+1,ndepth+1,4])
con = points_organized[:,:,:,-1]

##print('-------points_organized list--------')
##print(points_organized)

endtime1 = timeit.default_timer()
print("The time for data generation", endtime1 - starttime1)

starttime2 = timeit.default_timer()

connectivity = np.empty(shape=(ncells,8),dtype='uint32')
index_con = 0
for k in range(nrow):
    for j in range(ncolumn):
        for i in range(ndepth):
            arr1 = np.ravel(con[k,:,:][j:j+2,i:i+2], 'F')
            arr2 = np.ravel(con[k+1,:,:][j:j+2,i:i+2], 'F')
            vertex = np.array(list(zip(arr1,arr2)),dtype='uint32').flatten()          
            connectivity[index_con,:] = vertex
            index_con += 1

            if index_con % 6000 == 0:
                print(round(index_con/npoints*100,2),'%')

endtime2 = timeit.default_timer()
print("The time for for-loops", endtime2 - starttime2)

print('-------Vertex list--------')
print(connectivity)

starttime3 = timeit.default_timer()

# offsets: only for element type 11.
offsets = np.arange(8,ncells*8+0.5,8, dtype='uint32')

# cellTypes types

cellType = vtk.vtkVoxel().GetCellType() # # get cell type for voxel
cellTypes = np.zeros((nrow, ncolumn, ndepth), dtype='uint32').flatten() + cellType



# dataArray

pressure = np.random.rand(ncells).reshape( (nrow, ncolumn, ndepth))
temp = np.random.rand(npoints).reshape( (nrow + 1, ncolumn + 1, ndepth + 1))


##################################
### create xml-based vtu file  ###
##################################

PolyDMT = minidom.Document()

root = PolyDMT.createElement('VTKFile')
root.setAttribute('type', 'UnstructuredGrid')  # 属性是加在标签名之后、尖括号<>之中的
root.setAttribute('version', '0.1')
root.setAttribute('byte_order', 'LittleEndian')
PolyDMT.appendChild(root)

GridTypeChild = PolyDMT.createElement('UnstructuredGrid')
root.appendChild(GridTypeChild)

piece = PolyDMT.createElement('Piece')
piece.setAttribute('NumberOfPoints', str(npoints))
piece.setAttribute('NumberOfCells', str(ncells))
GridTypeChild.appendChild(piece)

Points = PolyDMT.createElement('Points')
piece.appendChild(Points)

cells = PolyDMT.createElement('Cells')
piece.appendChild(cells)

pointData = PolyDMT.createElement('PointData')
piece.appendChild(pointData)

cellData = PolyDMT.createElement('CellData')
piece.appendChild(cellData)

a = np.arange(0, 10)


def np2str(np_arr):
    np_arr = np_arr.flatten()
    str_arr = ""
    for i in np_arr:
        str_arr += str(i) + " "
    return str_arr


def dataFun(type, noComponents, dataContent, appendNode, name='', format='ascii'):
    '''
    创建元素、属性及其dataArray
    '''
    dataArray = PolyDMT.createElement('DataArray')
    if appendNode.tagName != "Points": dataArray.setAttribute('Name', name)
    dataArray.setAttribute('type', type)
    if appendNode.tagName != "Cells": dataArray.setAttribute('NumberOfComponents', noComponents)
    dataArray.setAttribute('format', format)
    appendNode.appendChild(dataArray)
    dataArray.appendChild(PolyDMT.createTextNode(np2str(dataContent)))


dataFun('Float32', '3', points[:,:3][:npoints], Points)
dataFun('Int32', '1', connectivity[:ncells*8], cells, 'connectivity')
dataFun('Int32', '1', offsets[:ncells], cells, 'offsets')
dataFun('UInt8', '1', cellTypes[:ncells], cells, 'types')

#greylevel = np.load('./greylevel.npy')

##dataFun('Float32', '1', greylevel, cellData, 'greylevel')
##dataFun('Int32', '3', a, cellData, 'YarnTangent')
##dataFun('Int32', '2', a, cellData, 'Location')
##dataFun('Int32', '1', a, cellData, 'VolumeFraction')
##dataFun('Int32', '1', a, cellData, 'SurfaceDistance')
dataFun('Float32', '1', temp, pointData, 'Temperature')

# write out the xml file.
xml_str = PolyDMT.toprettyxml(indent="\t")
save_path_file = "./createPolyDMT30.vtu"
with open(save_path_file, "w", buffering = 1) as f:
    f.write(xml_str)
f.close()
endtime3 = timeit.default_timer()

print("The time for file output", endtime3 - starttime3)

2. xml文件树的创建

import xml.dom.minidom as minidom
import numpy as np
import xml.sax

# class VtkHandler(xml.sax.ContentHandler):
#     def startElement(self, name, attrs):
#         print(name, attrs)
# handler = VtkHandler()
# parser = xml.sax.make_parser()
# parser.setContentHandler(handler)
# parser.parse('./texgen.vtu')

# domtree = xml.dom.minidom.parse('./texgen.vtu')
#
# group= domtree.documentElement
#
# for i in domtree.childNodes:
#     print(i.tagName)
# VTKFile = group.getElementsByTagName("VTKFile")
# print(VTKFile)
#
# DataArrays = group.getElementsByTagName('DataArray')
# # print("Name:{} ".format(group.getElementsByTagName('DataArray')[10].childNodes[0].data)) # 输出节点的数值
# # print("%d DataArrays" % DataArrays.length)
# for DataArray in DataArrays:
#     if DataArray.hasAttribute('Name'):
#         print('---'*200)
#         print(DataArray.getAttribute('Name'))
#         print(DataArray.childNodes[0].data)

##
### 修改原有标签的值
##persons[0].setAttribute("id","100")
##persons[1].getElementsByTagName('name')[0].childNodes[0].nodeValue = "New Name"
##persons[2].getElementsByTagName('age')[0].childNodes[0].nodeValue = "-10"
##


from xml.dom import minidom
import os

PolyDMT = minidom.Document()

root = PolyDMT.createElement('VTKFile')
root.setAttribute('type', 'UnstructuredGrid')  # 属性是加在标签名之后、尖括号<>之中的
root.setAttribute('version', '0.1')
root.setAttribute('byte_order', 'LittleEndian')
PolyDMT.appendChild(root)

GridTypeChild = PolyDMT.createElement('UnstructuredGrid')
root.appendChild(GridTypeChild)

piece = PolyDMT.createElement('Piece')
piece.setAttribute('NumberOfPoints', '216')
piece.setAttribute('NumberOfCells', '125')
GridTypeChild.appendChild(piece)

points = PolyDMT.createElement('Points')
piece.appendChild(points)

cells = PolyDMT.createElement('Cells')
piece.appendChild(cells)

pointData = PolyDMT.createElement('PointData')
piece.appendChild(pointData)

cellData = PolyDMT.createElement('CellData')
piece.appendChild(cellData)

a = np.arange(0, 10)


def np2str(np_arr):
    np_arr = np_arr.reshape(1, -1)[0]
    str_arr = ""
    for i in a:
        str_arr += str(i) + " "
    return str_arr

print(points.tagName=='Points')


def dataFun(type, noComponents, dataContent, appendNode, name='', format='ascii'):
    dataArray = PolyDMT.createElement('DataArray')
    if appendNode.tagName != "Points": dataArray.setAttribute('Name', name)
    dataArray.setAttribute('type', type)
    if appendNode.tagName != "Cells": dataArray.setAttribute('NumberOfComponents', noComponents)
    dataArray.setAttribute('format', format)
    appendNode.appendChild(dataArray)
    dataArray.appendChild(PolyDMT.createTextNode(np2str(dataContent)))

dataFun('Float32', '3', a, points)
dataFun('Int32', '1', a, cells, 'connectivity')
dataFun('Int32', '1', a, cells, 'offsets')
dataFun('Int32', '1', a, cells, 'types')
dataFun('Int32', '1', a,  cellData, 'YarnIndex')
dataFun('Int32', '3', a,  cellData, 'YarnTangent')
dataFun('Int32', '2', a,  cellData, 'Location')
dataFun('Int32', '1', a,  cellData, 'VolumeFraction')
dataFun('Int32', '1', a,  cellData, 'SurfaceDistance')
dataFun('Int32', '3', a,  cellData, 'Orientation')

# write out the xml file.
xml_str = PolyDMT.toprettyxml(indent="\t")
save_path_file = "data1.vtu"
with open(save_path_file, "w") as f:
    f.write(xml_str)

3. 将`TexGen`生成的体素网格转换为图片序列

这一功能已经完全重写并集成在了PolyKriging包中。引用包的函数的方法如下：

import pyvista as pv
import polykriging as pk
mesh = pv.read("./v2i.vtu")
mesh_shape = (20, 20, 5) # number of cells in x, y and z direction
# convert the mesh to images
pk.io.voxel2img(mesh, mesh_shape, dataset="YarnIndex", save_path="./img/", 
                scale=50, img_name="img", format="tif", 
                scale_algrithm="linear")

import xml.dom.minidom
import numpy as np
from skimage import io
from mayavi import mlab

# 将TexGen生成的模型转换为图片
# 注意像素尺寸

nx, ny, nz = 50, 50, 20  # number of cells in each direction

# load the unstructed grid (voxel model from TexGen with extension of .vtu)
domtree = xml.dom.minidom.parse('C:/Users/palme/Desktop/111.vtu')
rootNode = domtree.documentElement
UnstructuredGrid = rootNode.getElementsByTagName('UnstructuredGrid')

Piece = UnstructuredGrid[0].getElementsByTagName('Piece')
childlist = Piece[0].childNodes

CellData = Piece[0].getElementsByTagName('CellData')
DataArray = CellData[0].getElementsByTagName('DataArray')

###################################################
# Add grayscale as mask (field data assigned to each cell)
##io.use_plugin('pil') 
##ic = io.imread_collection(
##        'D:/04_coding/Python/06_Visualization/Mesh/00_data_acquisition/3D_2_4_5layers/Long_0/Long_0_img/*.tif') # 
####io.imshow_collection(ic)
####io.show()
##image_3d = io.concatenate_images(ic) # 

#3mlab.contour3d(image_3d[:,70:130,45:130])
#mlab.outline()
#np.save('./greylevel.npy',image_3d[:,70:130,45:130])


def np2str(np_arr):
        np_arr = np_arr.reshape(1,-1)[0]
        str_arr = ""
        for i in np_arr:
                str_arr += str(i)+" "
        return str_arr

#交换数组的轴至：shape(z,y,x)
image_3d = np.swapaxes(image_3d, 0,1)
image_3d = np.swapaxes(image_3d, 1,2)

imgId= domtree.createElement('DataArray')
# 切片时一定要注意多余数据不影响显示，故ParaView不会报错
imgId.appendChild(domtree.createTextNode(np2str(image_3d[:, :, 0:704:4])))
imgId.setAttribute('Name','image_3d')
imgId.setAttribute('NumberOfComponents','1')
imgId.setAttribute('format','ascii')
imgId.setAttribute('type','Float64')
CellData[0].appendChild(imgId)

# update DataArray
DataArray = CellData[0].getElementsByTagName('DataArray')
# Save vtu as image sequence
for item in DataArray:
        if item.getAttribute('Name') == "YarnIndex":
                YarnIndex = item.childNodes[0].data
                YarnIndex=YarnIndex.split(" ")
                img_sequence = np.empty([nz, nx, ny])
                for i in range(nz):
                        start = i*nx*ny
                        end = (i+1)*nx*ny
                        img = np.array(YarnIndex[start:end], dtype = float).reshape([nx,ny])
                        img_sequence[i,:,:]  = img 
                        io.imsave("C:/Users/palme/Desktop/img/" + str(i)+'.jpeg', img)

# Delete unnecessary cell data
for item in DataArray:
        if item.getAttribute('Name') in ["SurfaceDistance", "Location", "YarnTangent", "Orientation"]:
                CellData[0].removeChild(item)



# write out the xml file.
domtree.writexml(open("C:/Users/palme/Desktop/result.vtu","w"))
### Write the result to a vtu file
##with open("C:/Users/palme/Desktop/result.vtu","w") as fs:
##    fs.write(domtree.toxml())
##    fs.close()

4. 相关文件：

上一篇文章

[旧代码] 生成空间平面

下一篇文章

images to latex

查看评论 - NOTHING

Comments NOTHING

暂无评论

取消回复

Markdown Supported while Forbidden

你是我一生只会遇见一次的惊喜 ...

戳我呀 OωO 嘿嘿嘿ヾ(≧∇≦*)ゝ

bilibili~	(=・ω・=)	Tieba