Converting VTK structured grids to matrices
Tags: programming, research
I recently had to convert several structured grids in VTK format for further processing. Since I could not find any other code for doing this, here is my quick and dirty approach:
#!/usr/bin/env python3 from PIL import Image import numpy import os import re import sys reDimensions = re.compile( r'DIMENSIONS\s+(\d+)\s+(\d+)\s+(\d+)' ) reSpacing = re.compile( r'SPACING\s+(\d+\.\d+)\s+(\d+\.\d+)\s+(\d+\.\d+)' ) rePointData = re.compile( r'POINT_DATA\s+(\d+)' ) reTableData = re.compile( r'(\d+\.\d+)\s+(\d+\.\d+)\s+(\d+\.\d+)' ) for filename in sys.argv[1:]: with open( filename ) as f: width = 0 height = 0 depth = 0 data = None # Contains matrix data ii = 0 # Current data index for filling matrix # Stores only those coordinates (given as an x and y position) # that are nonzero. This permits to match adjacent time steps # in a very concise and sparse manner. nonZero = list() for line in f: if reDimensions.match(line): result = reDimensions.match(line) width = int( result.group(1) ) height = int( result.group(2) ) depth = int( result.group(3) ) assert width * height * depth > 0, "Invalid dimensions" data = numpy.zeros( (height,width) ) elif reSpacing.match(line): result = reSpacing.match(line) sx = result.group(1) sy = result.group(2) sz = result.group(3) # TODO: I am ignoring the spacing, hence I am only # able to load uniform rectilinear grids. assert sx == sy and sy == sz, "Invalid spacing" elif rePointData.match(line): result = rePointData.match(line) num = int( result.group(1) ) assert num == width * height * depth, "Inconsistent dimensions" elif reTableData.match(line): result = reTableData.match(line) d = (u,v,w) = ( float( result.group(1) ), float( result.group(2) ), float( result.group(3) ) ) for k in range(0,3): ix = ii % width iy = ii // width ii = ii + 1 if d[k] > 0: nonZero.append( (ix,iy) ) data[iy][ix] = d[k] # TODO: `data` is now a matrix corresponding to the dimensions of # the original grid.Do something with it... # # Store non-zero coordinates # outputCoordinates = "/tmp/"\ + os.path.splitext( os.path.basename(filename) )\ + ".txt" with open(outputCoordinates, "w") as g: for (x,y) in nonZero: print("%d\t%d" % (x,y), file=g) # # Store image # outputImage = "/tmp/"\ + os.path.splitext( os.path.basename(filename) )\ + ".png"
The code from above iterates over a list of VTK files (assumed to
be in ASCII format) that contains
POINT_DATA, i.e. a scalar value for
each grid cell. Provided that the grid follows uniform spacing, the code
is capable of parsing it and stores the
POINT_DATA values in a
matrix. Moreover, the list of non-zero coordinates of the matrix is
/tmp, along with a pixel-based image of grid. Both outputs
are highly useful for debugging and understanding the input.
Feel free to do more with this code—I am hereby releasing it into the public domain.