diff --git a/processing/post/addGrainID.py b/processing/post/addGrainID.py deleted file mode 100755 index 6493736d8..000000000 --- a/processing/post/addGrainID.py +++ /dev/null @@ -1,176 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: UTF-8 no BOM -*- - -import os,sys,copy -import numpy as np -import damask -from optparse import OptionParser -from scipy import spatial - -scriptName = os.path.splitext(os.path.basename(__file__))[0] -scriptID = ' '.join([scriptName,damask.version]) - - -parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """ -Add grain index based on similiarity of crystal lattice orientation. - -""", version = scriptID) - -parser.add_option('-r', - '--radius', - dest = 'radius', - type = 'float', metavar = 'float', - help = 'search radius') -parser.add_option('-d', - '--disorientation', - dest = 'disorientation', - type = 'float', metavar = 'float', - help = 'disorientation threshold in degrees [%default]') -parser.add_option('-s', - '--symmetry', - dest = 'symmetry', type = 'choice', choices = damask.Symmetry.lattices[1:], - metavar = 'string', - help = 'crystal symmetry [%default] {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:]))) -parser.add_option('-o', - '--orientation', - dest = 'quaternion', - metavar = 'string', - help = 'label of crystal orientation given as unit quaternion [%default]') -parser.add_option('-p', - '--pos', '--position', - dest = 'pos', - metavar = 'string', - help = 'label of coordinates [%default]') -parser.add_option('--quiet', - dest='verbose', - action = 'store_false', - help = 'hide status bar (useful when piping to file)') - -parser.set_defaults(disorientation = 5, - verbose = True, - quaternion = 'orientation', - symmetry = damask.Symmetry.lattices[-1], - pos = 'pos', - ) - -(options, filenames) = parser.parse_args() - -if options.radius is None: - parser.error('no radius specified.') - -cos_disorientation = np.cos(np.radians(options.disorientation/2.)) # cos of half the disorientation angle - -# --- loop over input files ------------------------------------------------------------------------- - -if filenames == []: filenames = [None] - -for name in filenames: - try: table = damask.ASCIItable(name = name, - buffered = False) - except: continue - damask.util.report(scriptName,name) - -# ------------------------------------------ read header ------------------------------------------- - - table.head_read() - -# ------------------------------------------ sanity checks ----------------------------------------- - - errors = [] - remarks = [] - - if not 3 >= table.label_dimension(options.pos) >= 1: - errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos)) - if not np.all(table.label_dimension(options.quaternion) == 4): - errors.append('input "{}" does not have dimension 4.'.format(options.quaternion)) - else: column = table.label_index(options.quaternion) - - if remarks != []: damask.util.croak(remarks) - if errors != []: - damask.util.croak(errors) - table.close(dismiss = True) - continue - -# ------------------------------------------ assemble header --------------------------------------- - - table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) - table.labels_append('grainID_{}@{:g}'.format(options.quaternion,options.disorientation)) # report orientation source and disorientation - table.head_write() - -# ------------------------------------------ build KD tree ----------------------------------------- - - table.data_readArray(options.pos) # read position vectors - grainID = -np.ones(len(table.data),dtype=int) - Npoints = table.data.shape[0] - kdtree = spatial.KDTree(copy.deepcopy(table.data)) - -# ------------------------------------------ assign grain IDs -------------------------------------- - - orientations = [] # quaternions found for grain - memberCounts = [] # number of voxels in grain - p = 0 # point counter - g = 0 # grain counter - matchedID = -1 - lastDistance = np.dot(kdtree.data[-1]-kdtree.data[0],kdtree.data[-1]-kdtree.data[0]) # (arbitrarily) use diagonal of cloud - - table.data_rewind() - while table.data_read(): # read next data line of ASCII table - - if options.verbose and Npoints > 100 and p%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero - damask.util.progressBar(iteration=p,total=Npoints) - - o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))), - symmetry = options.symmetry).reduced() - - matched = False - alreadyChecked = {} - candidates = [] - bestDisorientation = damask.Quaternion([0,0,0,1]) # initialize to 180 deg rotation as worst case - - for i in kdtree.query_ball_point(kdtree.data[p],options.radius): # check all neighboring points - gID = grainID[i] - if gID != -1 and gID not in alreadyChecked: # indexed point belonging to a grain not yet tested? - alreadyChecked[gID] = True # remember not to check again - disorientation = o.disorientation(orientations[gID],SST = False)[0] # compare against other orientation - if disorientation.quaternion.q > cos_disorientation: # within threshold ... - candidates.append(gID) # remember as potential candidate - if disorientation.quaternion.q >= bestDisorientation.q: # ... and better than current best? - matched = True - matchedID = gID # remember that grain - bestDisorientation = disorientation.quaternion - - if matched: # did match existing grain - memberCounts[matchedID] += 1 - if len(candidates) > 1: # ambiguity in grain identification? - largestGrain = sorted(candidates,key=lambda x:memberCounts[x])[-1] # find largest among potential candidate grains - matchedID = largestGrain - for c in [c for c in candidates if c != largestGrain]: # loop over smaller candidates - memberCounts[largestGrain] += memberCounts[c] # reassign member count of smaller to largest - memberCounts[c] = 0 - grainID = np.where(np.in1d(grainID,candidates), largestGrain, grainID) # relabel grid points of smaller candidates as largest one - - else: # no match -> new grain found - orientations += [o] # initialize with current orientation - memberCounts += [1] # start new membership counter - matchedID = g - g += 1 # increment grain counter - - grainID[p] = matchedID # remember grain index assigned to point - p += 1 # increment point - - grainIDs = np.where(np.array(memberCounts) > 0)[0] # identify "live" grain identifiers - packingMap = dict(zip(list(grainIDs),range(len(grainIDs)))) # map to condense into consecutive IDs - - table.data_rewind() - - outputAlive = True - p = 0 - damask.util.progressBar(iteration=1,total=1) - while outputAlive and table.data_read(): # read next data line of ASCII table - table.data_append(1+packingMap[grainID[p]]) # add (condensed) grain ID - outputAlive = table.data_write() # output processed line - p += 1 - -# ------------------------------------------ output finalization ----------------------------------- - - table.close() # close ASCII tables