474 lines
19 KiB
Python
Executable File
474 lines
19 KiB
Python
Executable File
#!/usr/bin/env python2.7
|
|
# -*- coding: UTF-8 no BOM -*-
|
|
|
|
import os,string,math,sys
|
|
import numpy as np
|
|
from optparse import OptionParser
|
|
import vtk
|
|
import damask
|
|
|
|
scriptName = os.path.splitext(os.path.basename(__file__))[0]
|
|
scriptID = ' '.join([scriptName,damask.version])
|
|
|
|
# -----------------------------
|
|
def getHeader(filename,sizeFastIndex,sizeSlowIndex,stepsize):
|
|
"""Returns header for ang file step size in micrometer"""
|
|
return '\n'.join([ \
|
|
'# TEM_PIXperUM 1.000000', \
|
|
'# x-star 1.000000', \
|
|
'# y-star 1.000000', \
|
|
'# z-star 1.000000', \
|
|
'# WorkingDistance 18.000000', \
|
|
'#', \
|
|
'# Phase 1', \
|
|
'# MaterialName XX', \
|
|
'# Formula XX', \
|
|
'# Info', \
|
|
'# Symmetry 43', \
|
|
'# LatticeConstants 2.870 2.870 2.870 90.000 90.000 90.000', \
|
|
'# NumberFamilies 1', \
|
|
'# hklFamilies 1 1 0 1 0.000000 1', \
|
|
'# Categories 0 0 0 0 0 ', \
|
|
'#', \
|
|
'# GRID: SqrGrid', \
|
|
'# XSTEP: ' + str(stepsize*1e6), \
|
|
'# YSTEP: ' + str(stepsize*1e6), \
|
|
'# NCOLS_ODD: ' + str(sizeFastIndex), \
|
|
'# NCOLS_EVEN: ' + str(sizeFastIndex), \
|
|
'# NROWS: ' + str(sizeSlowIndex), \
|
|
'#', \
|
|
'# OPERATOR: ' + string.replace('$Id$','\n','\\n'), \
|
|
'#', \
|
|
'# SAMPLEID: %s'%filename, \
|
|
'#', \
|
|
'# SCANID: ', \
|
|
'#', \
|
|
]) + '\n'
|
|
|
|
|
|
# -----------------------------
|
|
def positiveRadians(angle):
|
|
"""Returns positive angle in radians from angle in degrees"""
|
|
angle = math.radians(float(angle))
|
|
while angle < 0.0:
|
|
angle += 2.0 * math.pi
|
|
|
|
return angle
|
|
|
|
|
|
# -----------------------------
|
|
def getDataLine(angles,x,y,validData=True):
|
|
"""
|
|
Returns string of one line in ang file.
|
|
|
|
Convention in ang file: y coordinate comes first and is fastest index
|
|
positions in micrometer
|
|
"""
|
|
info = {True: (9999.9, 1.0, 0,99999,0.0),
|
|
False: ( -1.0,-1.0,-1, -1,1.0)}
|
|
return '%9.5f %9.5f %9.5f %12.5f %12.5f %6.1f %6.3f %2i %6i %6.3f \n'\
|
|
%(tuple(map(positiveRadians,angles))+(y*1e6,x*1e6)+info[validData])
|
|
|
|
|
|
|
|
# --------------------------------------------------------------------
|
|
# MAIN FUNCTION STARTS HERE
|
|
# --------------------------------------------------------------------
|
|
|
|
parser = OptionParser(usage='%prog options [file[s]]', description = """
|
|
Builds a ang files from a vtk file.
|
|
|
|
""", version = scriptID)
|
|
|
|
|
|
parser.add_option('--disp','--displacement',dest='dispLabel',
|
|
metavar ='string',
|
|
help='label of displacements [%default]')
|
|
parser.add_option('--euler', dest='eulerLabel', nargs=3,
|
|
metavar ='string string string',
|
|
help='labels of euler angles [%default]')
|
|
parser.add_option('-n','--normal', dest='normal', type='float', nargs=3,
|
|
metavar ='float float float',
|
|
help='normal of slices in direction of increasing slice numbers [%default]')
|
|
parser.add_option('-u','--up', dest='up', type='float', nargs=3,
|
|
metavar ='float float float',
|
|
help='up direction of slices [%default]')
|
|
parser.add_option('-i','--slices', dest='Nslices', type='int',
|
|
metavar ='int',
|
|
help='number of slices [%default]')
|
|
parser.add_option('-d','--distance', dest='distance', type='float',
|
|
metavar ='float',
|
|
help='slice distance [%default]')
|
|
parser.add_option('-s','--scale', dest='scale', type='float',
|
|
metavar ='float',
|
|
help='scale length from vtk file [%default]')
|
|
parser.add_option('-r','--resolution', dest='resolution', type='float',
|
|
metavar ='float',
|
|
help='scaling factor for resolution [%default]')
|
|
parser.add_option('--hex','--hexagonal', dest='hexagonal', action='store_true',
|
|
help='use in plane hexagonal grid')
|
|
parser.add_option('--interpolation', dest='interpolation', type='int',
|
|
metavar='float',
|
|
help='number of points for linear interpolation [%default]')
|
|
parser.add_option('--verbose', dest='verbose', action='store_true',
|
|
help='verbose mode')
|
|
parser.add_option('--visualize', dest='visualize', action='store_true',
|
|
help='visualize geometry')
|
|
|
|
parser.set_defaults(dispLabel = 'displacement')
|
|
parser.set_defaults(eulerLabel = ['1_1_eulerangles','1_2_eulerangles','1_3_eulerangles'])
|
|
parser.set_defaults(hexagonal = False)
|
|
parser.set_defaults(normal = [0.0,0.0,-1.0])
|
|
parser.set_defaults(up = [0.0,1.0,0.0])
|
|
parser.set_defaults(Nslices = 1)
|
|
parser.set_defaults(distance = 0.0)
|
|
parser.set_defaults(scale = 1.0)
|
|
parser.set_defaults(resolution = 1.0)
|
|
parser.set_defaults(dispScaling = 1.0)
|
|
parser.set_defaults(interpolation = 1)
|
|
parser.set_defaults(verbose = False)
|
|
parser.set_defaults(visualize = False)
|
|
(options,filenames) = parser.parse_args()
|
|
|
|
|
|
#--- SANITY CHECKS
|
|
|
|
# check for valid filenames
|
|
|
|
for filename in filenames:
|
|
if not os.path.exists(filename):
|
|
parser.error('file "%s" does not exist'%filename)
|
|
if not os.path.splitext(filename)[1] == '.vtk':
|
|
parser.error('"%s": need vtk file'%filename)
|
|
|
|
|
|
# check for othogonality of normal and up vector
|
|
|
|
if np.dot(np.array(options.normal),np.array(options.up)) > 1e-3:
|
|
parser.error('normal vector and up vector have to be orthogonal')
|
|
|
|
|
|
# check for options that are not yet implemented
|
|
|
|
if options.interpolation > 1:
|
|
parser.error('interpolation not yet supported')
|
|
if options.hexagonal:
|
|
parser.error('hexagonal grid not yet supported')
|
|
|
|
|
|
|
|
#--- ITERATE OVER FILES AND PROCESS THEM
|
|
|
|
for filename in filenames:
|
|
|
|
if options.verbose: sys.stdout.write("\nREADING VTK FILE\n")
|
|
# Read the source file
|
|
reader = vtk.vtkUnstructuredGridReader()
|
|
reader.SetFileName(filename)
|
|
reader.ReadAllScalarsOn()
|
|
reader.ReadAllVectorsOn()
|
|
reader.Update()
|
|
undeformedMesh = reader.GetOutput()
|
|
|
|
|
|
# Get euler angles from cell data
|
|
|
|
if options.verbose: sys.stdout.write("\nGETTING EULER ANGLES\n")
|
|
angles = {}
|
|
for i in range(reader.GetNumberOfScalarsInFile()):
|
|
scalarName = reader.GetScalarsNameInFile(i)
|
|
if scalarName in options.eulerLabel:
|
|
angles[scalarName] = undeformedMesh.GetCellData().GetScalars(scalarName)
|
|
if options.verbose: sys.stdout.write(" found scalar with name %s\n"%scalarName)
|
|
if len(angles) < 3: # found data for all three euler angles?
|
|
for label in options.eulerLabel:
|
|
if label not in angles.keys():
|
|
parser.error('Could not find scalar data with name %s'%label)
|
|
|
|
|
|
# Get deformed mesh
|
|
|
|
if options.verbose: sys.stdout.write("\nDEFORM MESH\n")
|
|
warpVector = vtk.vtkWarpVector()
|
|
undeformedMesh.GetPointData().SetActiveVectors(options.dispLabel)
|
|
warpVector.SetInput(undeformedMesh)
|
|
warpVector.Update()
|
|
deformedMesh = warpVector.GetOutput()
|
|
box = deformedMesh.GetBounds() # bounding box in mesh system
|
|
if options.verbose:
|
|
sys.stdout.write(" bounding box in lab system\n")
|
|
sys.stdout.write(" x (% .8f % .8f)\n"%(box[0],box[1]))
|
|
sys.stdout.write(" y (% .8f % .8f)\n"%(box[2],box[3]))
|
|
sys.stdout.write(" z (% .8f % .8f)\n"%(box[4],box[5]))
|
|
|
|
|
|
# Get cell centers of deformed mesh (position of ips)
|
|
|
|
if options.verbose: sys.stdout.write("\nGETTING CELL CENTERS OF DEFORMED MESH\n")
|
|
cellCenter = vtk.vtkCellCenters()
|
|
cellCenter.SetVertexCells(0) # do not generate vertex cells, just points
|
|
cellCenter.SetInput(deformedMesh)
|
|
cellCenter.Update()
|
|
meshIPs = cellCenter.GetOutput()
|
|
|
|
|
|
# Get outer surface of deformed mesh
|
|
|
|
if options.verbose: sys.stdout.write("\nGETTING OUTER SURFACE OF DEFORMED MESH\n")
|
|
surfaceFilter = vtk.vtkDataSetSurfaceFilter()
|
|
surfaceFilter.SetInput(deformedMesh)
|
|
surfaceFilter.Update()
|
|
surface = surfaceFilter.GetOutput()
|
|
|
|
|
|
# Get coordinate system for ang files
|
|
# z-vector is normal to slices
|
|
# x-vector corresponds to the up-direction
|
|
# "R" rotates coordinates from the mesh system into the TSL system
|
|
|
|
if options.verbose: sys.stdout.write("\nGETTING COORDINATE SYSTEM FOR ANG FILES\n")
|
|
z = np.array(options.normal,dtype='float')
|
|
z = z / np.linalg.norm(z)
|
|
x = np.array(options.up,dtype='float')
|
|
x = x / np.linalg.norm(x)
|
|
y = np.cross(z,x)
|
|
R = np.array([x,y,z])
|
|
if options.verbose:
|
|
sys.stdout.write(" axis (x: up direction, z: slice normal)\n")
|
|
sys.stdout.write(" x (% .8f % .8f % .8f)\n"%tuple(x))
|
|
sys.stdout.write(" y (% .8f % .8f % .8f)\n"%tuple(y))
|
|
sys.stdout.write(" z (% .8f % .8f % .8f)\n"%tuple(z))
|
|
|
|
|
|
# Get bounding box in rotated system (x,y,z)
|
|
|
|
if options.verbose: sys.stdout.write("\nGETTING BOUNDING BOX IN ROTATED SYSTEM\n")
|
|
rotatedbox = [[np.inf,-np.inf] for i in range(3)] # bounding box in rotated TSL system
|
|
for n in range(8): # loop over eight vertices of mesh bounding box
|
|
vert = np.array([box[0+(n/1)%2],
|
|
box[2+(n/2)%2],
|
|
box[4+(n/4)%2]]) # vertex in mesh system
|
|
rotatedvert = np.dot(R,vert) # vertex in rotated system
|
|
for i in range(3):
|
|
rotatedbox[i][0] = min(rotatedbox[i][0],rotatedvert[i])
|
|
rotatedbox[i][1] = max(rotatedbox[i][1],rotatedvert[i])
|
|
if options.verbose:
|
|
sys.stdout.write(" bounding box in rotated system\n")
|
|
sys.stdout.write(" x (% .8f % .8f)\n"%tuple(rotatedbox[0]))
|
|
sys.stdout.write(" y (% .8f % .8f)\n"%tuple(rotatedbox[1]))
|
|
sys.stdout.write(" z (% .8f % .8f)\n"%tuple(rotatedbox[2]))
|
|
|
|
|
|
# Correct bounding box so that a multiplicity of the resolution fits into it
|
|
# and get number of points and extent in each (rotated) axis direction
|
|
|
|
if options.verbose: sys.stdout.write("\nCORRECTING EXTENT OF BOUNDING BOX IN ROTATED SYSTEM\n")
|
|
correction = []
|
|
Npoints = []
|
|
extent = [rotatedbox[i][1] - rotatedbox[i][0] for i in range(3)]
|
|
for i in range(2):
|
|
Npoints.extend([int(math.ceil(extent[i] / options.resolution))])
|
|
correction.extend([float(Npoints[i]) * options.resolution - extent[i]])
|
|
if options.distance > 0.0:
|
|
Npoints.extend([int(math.ceil(extent[2] / options.distance))])
|
|
correction.extend([float(Npoints[2]) * options.distance - extent[2]])
|
|
else:
|
|
Npoints.extend([options.Nslices])
|
|
correction.extend([0.0])
|
|
options.distance = extent[2] / float(options.Nslices)
|
|
for i in range(3):
|
|
rotatedbox[i][0] = rotatedbox[i][0] - 0.5 * correction[i]
|
|
rotatedbox[i][1] = rotatedbox[i][1] + 0.5 * correction[i]
|
|
extent[i] = rotatedbox[i][1] - rotatedbox[i][0]
|
|
NpointsPerSlice = Npoints[0] * Npoints[1]
|
|
totalNpoints = NpointsPerSlice * Npoints[2]
|
|
if options.verbose:
|
|
sys.stdout.write(" corrected bounding box in rotated system\n")
|
|
sys.stdout.write(" x (% .8f % .8f)\n"%tuple(rotatedbox[0]))
|
|
sys.stdout.write(" y (% .8f % .8f)\n"%tuple(rotatedbox[1]))
|
|
sys.stdout.write(" z (% .8f % .8f)\n"%tuple(rotatedbox[2]))
|
|
|
|
|
|
# Generate new regular point grid for ang files
|
|
# Use "polydata" object with points as single vertices
|
|
# beware of TSL convention: y direction is fastest index
|
|
|
|
if options.verbose: sys.stdout.write("\nGENERATING POINTS FOR POINT GRID")
|
|
points = vtk.vtkPoints()
|
|
for k in range(Npoints[2]):
|
|
for j in range(Npoints[0]):
|
|
for i in range(Npoints[1]): # y is fastest index
|
|
rotatedpoint = np.array([rotatedbox[0][0] + (float(j) + 0.5) * options.resolution,
|
|
rotatedbox[1][0] + (float(i) + 0.5) * options.resolution,
|
|
rotatedbox[2][0] + (float(k) + 0.5) * options.distance ]) # point in rotated system
|
|
point = np.dot(R.T,rotatedpoint) # point in mesh system
|
|
points.InsertNextPoint(list(point))
|
|
if options.verbose:
|
|
sys.stdout.write("\rGENERATING POINTS FOR POINT GRID %d%%" %(100*(Npoints[1]*(k*Npoints[0]+j)+i+1)/totalNpoints))
|
|
sys.stdout.flush()
|
|
if options.verbose:
|
|
sys.stdout.write("\n number of slices: %i\n"%Npoints[2])
|
|
sys.stdout.write(" slice spacing: %.8f\n"%options.distance)
|
|
if Npoints[2] > 1:
|
|
sys.stdout.write(" number of points per slice: %i = %i rows * %i points in row\n"%(NpointsPerSlice,Npoints[0],Npoints[1]))
|
|
sys.stdout.write(" grid resolution: %.8f\n"%options.resolution)
|
|
|
|
if options.verbose: sys.stdout.write("\nGENERATING VERTICES FOR POINT GRID")
|
|
vertices = vtk.vtkCellArray()
|
|
for i in range(totalNpoints):
|
|
vertex = vtk.vtkVertex()
|
|
vertex.GetPointIds().SetId(0,i) # each vertex consists of exactly one (index 0) point with ID "i"
|
|
vertices.InsertNextCell(vertex)
|
|
if options.verbose:
|
|
sys.stdout.write("\rGENERATING VERTICES FOR POINT GRID %d%%" %(100*(i+1)/totalNpoints))
|
|
sys.stdout.flush()
|
|
|
|
if options.verbose: sys.stdout.write("\n\nGENERATING POINT GRID\n")
|
|
pointgrid = vtk.vtkPolyData()
|
|
pointgrid.SetPoints(points)
|
|
pointgrid.SetVerts(vertices)
|
|
pointgrid.Update()
|
|
|
|
|
|
# Find out which points reside inside mesh geometry
|
|
|
|
if options.verbose: sys.stdout.write("\nIDENTIFYING POINTS INSIDE MESH GEOMETRY\n")
|
|
enclosedPoints = vtk.vtkSelectEnclosedPoints()
|
|
enclosedPoints.SetSurface(surface)
|
|
enclosedPoints.SetInput(pointgrid)
|
|
enclosedPoints.Update()
|
|
|
|
|
|
# Build kdtree from mesh IPs and match mesh IPs to point grid
|
|
|
|
if options.verbose: sys.stdout.write("\nBUILDING MAPPING OF GRID POINTS")
|
|
kdTree = vtk.vtkKdTree()
|
|
kdTree.BuildLocatorFromPoints(meshIPs.GetPoints())
|
|
gridToMesh = []
|
|
ids = vtk.vtkIdList()
|
|
NenclosedPoints = 0
|
|
for i in range(pointgrid.GetNumberOfPoints()):
|
|
gridToMesh.append([])
|
|
if enclosedPoints.IsInside(i):
|
|
NenclosedPoints += 1
|
|
# here one could use faster(?) "FindClosestPoint" if only first nearest neighbor required
|
|
kdTree.FindClosestNPoints(options.interpolation,pointgrid.GetPoint(i),ids)
|
|
for j in range(ids.GetNumberOfIds()):
|
|
gridToMesh[-1].extend([ids.GetId(j)])
|
|
if options.verbose:
|
|
sys.stdout.write("\rBUILDING MAPPING OF GRID POINTS %d%%" %(100*(i+1)/totalNpoints))
|
|
sys.stdout.flush()
|
|
if options.verbose:
|
|
sys.stdout.write("\n Number of points inside mesh geometry %i\n"%NenclosedPoints)
|
|
sys.stdout.write(" Number of points outside mesh geometry %i\n"%(totalNpoints - NenclosedPoints))
|
|
|
|
|
|
|
|
# ITERATE OVER SLICES AND CREATE ANG FILE
|
|
|
|
if options.verbose:
|
|
sys.stdout.write("\nWRITING OUT ANG FILES\n")
|
|
sys.stdout.write(" scaling all length with %f\n"%options.scale)
|
|
x0,y0,z0 = np.dot(R,pointgrid.GetPoint(0)) # first point on slice defines origin
|
|
for sliceN in range(Npoints[2]):
|
|
|
|
# Open file and write header
|
|
|
|
angfilename = eval('"'+eval("'%%s_slice%%0%ii.ang'%(math.log10(Npoints[2])+1)")+'"%(os.path.splitext(filename)[0],sliceN+1)')
|
|
with open(angfilename,'w') as angfile:
|
|
if options.verbose: sys.stdout.write(" %s\n"%angfilename)
|
|
angfile.write(getHeader(filename,Npoints[1],Npoints[0],options.resolution*options.scale))
|
|
for i in range(sliceN*NpointsPerSlice,(sliceN+1)*NpointsPerSlice): # Iterate over points on slice
|
|
|
|
|
|
# Get euler angles of closest IDs
|
|
|
|
if enclosedPoints.IsInside(i):
|
|
phi = []
|
|
for j in range(len(gridToMesh[i])):
|
|
IP = gridToMesh[i][j]
|
|
phi.append([])
|
|
for k in range(3):
|
|
phi[-1].extend([angles[options.eulerLabel[k]].GetValue(IP)])
|
|
else:
|
|
phi = [[720,720,720]] # fake angles
|
|
|
|
|
|
# Interpolate Euler angle
|
|
# NOT YET IMPLEMENTED, simply take the nearest neighbors values
|
|
|
|
interpolatedPhi = phi[0]
|
|
|
|
|
|
# write data to ang file
|
|
|
|
x,y,z = np.dot(R,pointgrid.GetPoint(i)) # point in rotated TSL system
|
|
x -= x0 # first point on slice defines origin
|
|
y -= y0 # first point on slice defines origin
|
|
x *= options.scale
|
|
y *= options.scale
|
|
angfile.write(getDataLine(interpolatedPhi,x,y,enclosedPoints.IsInside(i)))
|
|
|
|
|
|
# Visualize slices
|
|
|
|
if options.visualize:
|
|
meshMapper = vtk.vtkDataSetMapper()
|
|
meshMapper.SetInput(surface)
|
|
meshMapper.ScalarVisibilityOff() # do not use scalar data for coloring
|
|
meshActor = vtk.vtkActor()
|
|
meshActor.SetMapper(meshMapper)
|
|
meshActor.GetProperty().SetOpacity(0.2)
|
|
meshActor.GetProperty().SetColor(1.0,1.0,0)
|
|
meshActor.GetProperty().BackfaceCullingOn()
|
|
# meshActor.GetProperty().SetEdgeColor(1,1,0.5)
|
|
# meshActor.GetProperty().EdgeVisibilityOn()
|
|
|
|
boxpoints = vtk.vtkPoints()
|
|
for n in range(8):
|
|
P = [rotatedbox[0][(n/1)%2],
|
|
rotatedbox[1][(n/2)%2],
|
|
rotatedbox[2][(n/4)%2]]
|
|
boxpoints.InsertNextPoint(list(np.dot(R.T,np.array(P))))
|
|
box = vtk.vtkHexahedron()
|
|
for n,i in enumerate([0,1,3,2,4,5,7,6]):
|
|
box.GetPointIds().SetId(n,i)
|
|
boxgrid = vtk.vtkUnstructuredGrid()
|
|
boxgrid.SetPoints(boxpoints)
|
|
boxgrid.InsertNextCell(box.GetCellType(), box.GetPointIds())
|
|
boxsurfaceFilter = vtk.vtkDataSetSurfaceFilter()
|
|
boxsurfaceFilter.SetInput(boxgrid)
|
|
boxsurfaceFilter.Update()
|
|
boxsurface = boxsurfaceFilter.GetOutput()
|
|
|
|
boxMapper = vtk.vtkDataSetMapper()
|
|
boxMapper.SetInput(boxsurface)
|
|
boxActor = vtk.vtkActor()
|
|
boxActor.SetMapper(boxMapper)
|
|
boxActor.GetProperty().SetLineWidth(2.0)
|
|
boxActor.GetProperty().SetRepresentationToWireframe()
|
|
|
|
gridMapper = vtk.vtkDataSetMapper()
|
|
gridMapper.SetInput(pointgrid)
|
|
gridActor = vtk.vtkActor()
|
|
gridActor.SetMapper(gridMapper)
|
|
gridActor.GetProperty().SetColor(0,0,0)
|
|
gridActor.GetProperty().SetPointSize(3)
|
|
|
|
|
|
renderer = vtk.vtkRenderer()
|
|
renderWindow = vtk.vtkRenderWindow()
|
|
renderWindow.FullScreenOn()
|
|
renderWindow.AddRenderer(renderer)
|
|
renderWindowInteractor = vtk.vtkRenderWindowInteractor()
|
|
renderWindowInteractor.SetRenderWindow(renderWindow)
|
|
renderer.AddActor(meshActor)
|
|
renderer.AddActor(boxActor)
|
|
renderer.AddActor(gridActor)
|
|
renderer.SetBackground(1,1,1)
|
|
|
|
renderWindow.Render()
|
|
renderWindowInteractor.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera())
|
|
renderWindowInteractor.Start()
|
|
|