From 0e93d51fed46a8a010ed5f21729fb32e388192f7 Mon Sep 17 00:00:00 2001
From: Pratheek Shanthraj
Date: Fri, 18 Jan 2013 11:42:27 +0000
Subject: [PATCH] adopted same algorithm as in addEuclideanDistance
---
processing/pre/geom_euclideanDistance.py | 340 ++++++++++++++---------
1 file changed, 215 insertions(+), 125 deletions(-)
diff --git a/processing/pre/geom_euclideanDistance.py b/processing/pre/geom_euclideanDistance.py
index 0fed81883..f8250f235 100755
--- a/processing/pre/geom_euclideanDistance.py
+++ b/processing/pre/geom_euclideanDistance.py
@@ -1,96 +1,52 @@
#!/usr/bin/env python
-# -*- coding: UTF-8 no BOM -*-
-import os,sys,string,re,numpy,skfmm
-from optparse import OptionParser, OptionGroup, Option, SUPPRESS_HELP
+import os,re,sys,math,numpy,string,damask
from scipy import ndimage
+from optparse import OptionParser, Option
# -----------------------------
-class extendedOption(Option):
+class extendableOption(Option):
# -----------------------------
# used for definition of new option parser action 'extend', which enables to take multiple option arguments
# taken from online tutorial http://docs.python.org/library/optparse.html
-
- ACTIONS = Option.ACTIONS + ("extend",)
- STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",)
- TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",)
- ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",)
+
+ ACTIONS = Option.ACTIONS + ("extend",)
+ STORE_ACTIONS = Option.STORE_ACTIONS + ("extend",)
+ TYPED_ACTIONS = Option.TYPED_ACTIONS + ("extend",)
+ ALWAYS_TYPED_ACTIONS = Option.ALWAYS_TYPED_ACTIONS + ("extend",)
- def take_action(self, action, dest, opt, value, values, parser):
- if action == "extend":
- lvalue = value.split(",")
- values.ensure_value(dest, []).extend(lvalue)
- else:
- Option.take_action(self, action, dest, opt, value, values, parser)
-
-
-
-def outStdout(cmd,locals):
- if cmd[0:3] == '(!)':
- exec(cmd[3:])
- elif cmd[0:3] == '(?)':
- cmd = eval(cmd[3:])
- print cmd
- else:
- print cmd
- return
-
-def outFile(cmd,locals):
- if cmd[0:3] == '(!)':
- exec(cmd[3:])
- elif cmd[0:3] == '(?)':
- cmd = eval(cmd[3:])
- locals['filepointer'].write(cmd+'\n')
- else:
- locals['filepointer'].write(cmd+'\n')
- return
-
-
-def output(cmds,locals,dest):
- for cmd in cmds:
- if isinstance(cmd,list):
- output(cmd,locals,dest)
+ def take_action(self, action, dest, opt, value, values, parser):
+ if action == "extend":
+ lvalue = value.split(",")
+ values.ensure_value(dest, []).extend(lvalue)
else:
- {\
- 'File': outFile,\
- 'Stdout': outStdout,\
- }[dest](str(cmd),locals)
- return
+ Option.take_action(self, action, dest, opt, value, values, parser)
-# +++++++++++++++++++++++++++++++++++++++++++++++++++
-def vtk_writeASCII_mesh(dim,res,origin,data):
-# +++++++++++++++++++++++++++++++++++++++++++++++++++
- """ function writes data array defined on a rectilinear grid """
- N = res[0]*res[1]*res[2]
-
- cmds = [\
- '# vtk DataFile Version 3.1',
- string.replace('powered by $Id: spectral_geomCheck.py 1575 2012-06-26 18:07:38Z MPIE\p.eisenlohr $','\n','\\n'),
- 'ASCII',
- 'DATASET RECTILINEAR_GRID',
- 'DIMENSIONS %i %i %i'%(res[0]+1,res[1]+1,res[2]+1),
- 'X_COORDINATES %i float'%(res[0]+1),
- ' '.join(map(str,[i*dim[0]/res[0]+origin[0] for i in range(res[0]+1)])),
- 'Y_COORDINATES %i float'%(res[1]+1),
- ' '.join(map(str,[i*dim[1]/res[1]+origin[1] for i in range(res[1]+1)])),
- 'Z_COORDINATES %i float'%(res[2]+1),
- ' '.join(map(str,[i*dim[2]/res[2]+origin[2] for i in range(res[2]+1)])),
- 'CELL_DATA %i'%N,
- ]
-
- for datatype in data:
- for item in data[datatype]:
- cmds += [\
- '%s %s float'%(datatype.upper()+{True:'',False:'S'}[datatype.lower().endswith('s')],item),
- 'LOOKUP_TABLE default',
- [[['\t'.join(map(str,data[datatype][item][:,j,k]))] for j in range(res[1])] for k in range(res[2])]
- ]
+def periodic_3Dpad(array, rimdim=(1,1,1)):
- return cmds
+ rimdim = numpy.array(rimdim,'i')
+ size = numpy.array(array.shape,'i')
+ padded = numpy.empty(size+2*rimdim,array.dtype)
+ padded[rimdim[0]:rimdim[0]+size[0],
+ rimdim[1]:rimdim[1]+size[1],
+ rimdim[2]:rimdim[2]+size[2]] = array
+ p = numpy.zeros(3,'i')
+ for side in xrange(3):
+ for p[(side+2)%3] in xrange(padded.shape[(side+2)%3]):
+ for p[(side+1)%3] in xrange(padded.shape[(side+1)%3]):
+ for p[side%3] in xrange(rimdim[side%3]):
+ spot = (p-rimdim)%size
+ padded[p[0],p[1],p[2]] = array[spot[0],spot[1],spot[2]]
+ for p[side%3] in xrange(rimdim[side%3]+size[side%3],size[side%3]+2*rimdim[side%3]):
+ spot = (p-rimdim)%size
+ padded[p[0],p[1],p[2]] = array[spot[0],spot[1],spot[2]]
+ return padded
-# ----------------------- MAIN -------------------------------
+# --------------------------------------------------------------------
+# MAIN
+# --------------------------------------------------------------------
identifiers = {
'resolution': ['a','b','c'],
@@ -103,28 +59,108 @@ mappings = {
'origin': lambda x: float(x),
}
-parser = OptionParser(option_class=extendedOption, usage='%prog [geomfile[s]]', description = """
-Produce Euclidean distance map from geom description
+features = [ \
+ {'aliens': 1, 'names': ['boundary','biplane'],},
+ {'aliens': 2, 'names': ['tripleline',],},
+ {'aliens': 3, 'names': ['quadruplepoint',],}
+ ]
-""" + string.replace('$Id: spectral_geomCheck.py 1575 2012-06-26 18:07:38Z MPIE\p.eisenlohr $','\n','\\n')
+neighborhoods = {
+ 'neumann':numpy.array([
+ [-1, 0, 0],
+ [ 1, 0, 0],
+ [ 0,-1, 0],
+ [ 0, 1, 0],
+ [ 0, 0,-1],
+ [ 0, 0, 1],
+ ]),
+ 'moore':numpy.array([
+ [-1,-1,-1],
+ [ 0,-1,-1],
+ [ 1,-1,-1],
+ [-1, 0,-1],
+ [ 0, 0,-1],
+ [ 1, 0,-1],
+ [-1, 1,-1],
+ [ 0, 1,-1],
+ [ 1, 1,-1],
+ [-1,-1, 0],
+ [ 0,-1, 0],
+ [ 1,-1, 0],
+ [-1, 0, 0],
+#
+ [ 1, 0, 0],
+ [-1, 1, 0],
+ [ 0, 1, 0],
+ [ 1, 1, 0],
+ [-1,-1, 1],
+ [ 0,-1, 1],
+ [ 1,-1, 1],
+ [-1, 0, 1],
+ [ 0, 0, 1],
+ [ 1, 0, 1],
+ [-1, 1, 1],
+ [ 0, 1, 1],
+ [ 1, 1, 1],
+ ])
+ }
+
+parser = OptionParser(option_class=extendableOption, usage='%prog options [file[s]]', description = """
+Produce geom files containing Euclidean distance to grain structural features:
+boundaries, triple lines, and quadruple points.
+
+""" + string.replace('$Id: addEuclideanDistance.py 2039 2012-12-19 14:50:45Z MPIE\p.shanthraj $','\n','\\n')
)
-(options, filenames) = parser.parse_args()
+parser.add_option('-t','--type', dest='type', action='extend', type='string', \
+ help='feature type (%s)'%(', '.join(map(lambda x:', '.join(x['names']),features))))
+parser.add_option('-n','--neighborhood', dest='neigborhood', action='store', type='string', \
+ help='type of neighborhood (%s)'%(', '.join(neighborhoods.keys())), \
+ metavar='')
+parser.add_option('-2', '--twodimensional', dest='twoD', action='store_true', \
+ help='output geom file with two-dimensional data arrangement')
+parser.set_defaults(type = [])
+parser.set_defaults(neighborhood = 'neumann')
+parser.set_defaults(twoD = False)
+
+(options,filenames) = parser.parse_args()
+
+options.neighborhood = options.neighborhood.lower()
+if options.neighborhood not in neighborhoods:
+ parser.error('unknown neighborhood %s!'%options.neighborhood)
+
+feature_list = []
+for i,feature in enumerate(features):
+ for name in feature['names']:
+ for type in options.type:
+ if name.startswith(type):
+ feature_list.append(i) # remember valid features
+ break
+
+print feature_list
# ------------------------------------------ setup file handles ---------------------------------------
files = []
if filenames == []:
- files.append({'name':'STDIN', 'input':sys.stdin})
+ files.append({'name':'STDIN',
+ 'input':sys.stdin,
+ 'output':sys.stdout,
+ 'croak':sys.stderr,
+ })
else:
for name in filenames:
if os.path.exists(name):
- files.append({'name':name, 'input':open(name)})
+ files.append({'name':name,
+ 'input':open(name),
+ 'output':[open(features[feature]['names'][0]+'_'+name,'w') for feature in feature_list],
+ 'croak':sys.stdout,
+ })
# ------------------------------------------ loop over input files ---------------------------------------
for file in files:
- if file['name'] != 'STDIN': print file['name']
+ if file['name'] != 'STDIN': file['croak'].write(file['name']+'\n')
# get labels by either read the first row, or - if keyword header is present - the last line of the header
@@ -132,7 +168,7 @@ for file in files:
m = re.search('(\d+)\s*head', firstline.lower())
if m:
headerlines = int(m.group(1))
- headers = [file['input'].readline() for i in range(headerlines)]
+ headers = [firstline]+[file['input'].readline() for i in range(headerlines)]
else:
headerlines = 1
headers = firstline
@@ -140,56 +176,110 @@ for file in files:
content = file['input'].readlines()
file['input'].close()
- info = {'resolution': [0,0,0],
- 'dimension': [0.0,0.0,0.0],
- 'origin': [0.0,0.0,0.0],
+ info = {'resolution': numpy.array([0,0,0]),
+ 'dimension': numpy.array([0.0,0.0,0.0]),
+ 'origin': numpy.array([0.0,0.0,0.0]),
+ 'homogenization': 1,
}
+
+ new_header = []
for header in headers:
headitems = map(str.lower,header.split())
- if headitems[0] in identifiers.keys():
- for i in xrange(len(identifiers[headitems[0]])):
- info[headitems[0]][i] = mappings[headitems[0]](headitems[headitems.index(identifiers[headitems[0]][i])+1])
+ if headitems[0] in mappings.keys():
+ if headitems[0] in identifiers.keys():
+ for i in xrange(len(identifiers[headitems[0]])):
+ info[headitems[0]][i] = \
+ mappings[headitems[0]](headitems[headitems.index(identifiers[headitems[0]][i])+1])
+ else:
+ info[headitems[0]] = mappings[headitems[0]](headitems[1])
- if info['resolution'] == [0,0,0]:
- print 'no resolution info found.'
- sys.exit(1)
- if info['dimension'] == [0.0,0.0,0.0]:
- print 'no dimension info found.'
- sys.exit(1)
+ if numpy.all(info['resolution'] == 0):
+ file['croak'].write('no resolution info found.\n')
+ continue
+ if numpy.all(info['dimension'] == 0.0):
+ file['croak'].write('no dimension info found.\n')
+ continue
- print 'resolution: %s'%(' x '.join(map(str,info['resolution'])))
- print 'dimension: %s'%(' x '.join(map(str,info['dimension'])))
- print 'origin: %s'%(' : '.join(map(str,info['origin'])))
+ file['croak'].write('resolution: %s\n'%(' x '.join(map(str,info['resolution']))) + \
+ 'dimension: %s\n'%(' x '.join(map(str,info['dimension']))) + \
+ 'origin: %s\n'%(' : '.join(map(str,info['origin']))) + \
+ 'homogenization: %i\n'%info['homogenization'])
- dx = info['dimension'][0]/info['resolution'][0]
-
- data = {'scalar':{'perimeter':numpy.zeros(info['resolution'],'i'),
- 'distance':numpy.zeros(info['resolution'],'i')}}
+ new_header.append("resolution\ta %i\tb %i\tc %i\n"%(
+ info['resolution'][0],
+ info['resolution'][1],
+ info['resolution'][2],))
+ new_header.append("dimension\tx %f\ty %f\tz %f\n"%(
+ info['dimension'][0],
+ info['dimension'][1],
+ info['dimension'][2],))
+ new_header.append("origin\tx %f\ty %f\tz %f\n"%(
+ info['origin'][0],
+ info['origin'][1],
+ info['origin'][2],))
+ new_header.append("homogenization\t%i\n"%info['homogenization'])
+
+ structure = numpy.zeros(info['resolution'],'i')
i = 0
for line in content:
for item in map(int,line.split()):
- data['scalar']['perimeter'][i%info['resolution'][0],(i/info['resolution'][0])%info['resolution'][1],i/info['resolution'][0]/info['resolution'][1]] = item
+ structure[i%info['resolution'][0],
+ (i/info['resolution'][0])%info['resolution'][1],
+ i/info['resolution'][0] /info['resolution'][1]] = item
i += 1
+
+ neighborhood = neighborhoods[options.neighborhood]
+ convoluted = numpy.empty([len(neighborhood)]+list(info['resolution']+2),'i')
+ microstructure = periodic_3Dpad(structure)
+
+ for i,p in enumerate(neighborhood):
+ stencil = numpy.zeros((3,3,3),'i')
+ stencil[1,1,1] = -1
+ stencil[p[0]+1,
+ p[1]+1,
+ p[2]+1] = 1
-# data['scalar']['perimeter'] = numpy.where(ndimage.morphology.grey_dilation(data['scalar']['perimeter'],size=(3,3,3))-data['scalar']['perimeter']>0,0,1)
+ convoluted[i,:,:,:] = ndimage.convolve(microstructure,stencil)
- FDstencil_x = numpy.zeros([3,3,3])
- FDstencil_x[:,1,1] = [-1,0,1]
- FDstencil_y = numpy.zeros([3,3,3])
- FDstencil_y[1,:,1] = [-1,0,1]
- FDstencil_z = numpy.zeros([3,3,3])
- FDstencil_z[1,1,:] = [-1,0,1]
- data['scalar']['perimeter'] = numpy.where(numpy.abs(ndimage.convolve(data['scalar']['perimeter'], FDstencil_x)) + numpy.abs(ndimage.convolve(data['scalar']['perimeter'], FDstencil_y)) + numpy.abs(ndimage.convolve(data['scalar']['perimeter'], FDstencil_z))>0,0,1)
- data['scalar']['distance'] = skfmm.distance(data['scalar']['perimeter'], dx=dx)
+ distance = numpy.ones((len(feature_list),info['resolution'][0],info['resolution'][1],info['resolution'][2]),'d')
- out = {}
- out['mesh'] = vtk_writeASCII_mesh(info['dimension'],info['resolution'],info['origin'],data)
+ convoluted = numpy.sort(convoluted,axis=0)
+ uniques = numpy.zeros(info['resolution'])
+ check = numpy.empty(info['resolution'])
+ check[:,:,:] = numpy.nan
+ for i in xrange(len(neighborhood)):
+ uniques += numpy.where(convoluted[i,1:-1,1:-1,1:-1] == check,0,1)
+ check = convoluted[i,1:-1,1:-1,1:-1]
+ for i,feature_id in enumerate(feature_list):
+ distance[i,:,:,:] = numpy.where(uniques > features[feature_id]['aliens'],0.0,1.0)
- for what in out.keys():
- if file['name'] == 'STDIN':
- output(out[what],{},'Stdout')
- else:
- (head,tail) = os.path.split(file['name'])
- vtk = open(os.path.join(head,what+'_'+os.path.splitext(tail)[0]+'.vtk'), 'w')
- output(out[what],{'filepointer':vtk},'File')
- vtk.close()
+ for i in xrange(len(feature_list)):
+ distance[i,:,:,:] = ndimage.morphology.distance_transform_edt(distance[i,:,:,:])*[max(info['dimension']/info['resolution'])]*3
+
+
+ for i,feature in enumerate(feature_list):
+ formatwidth = int(math.floor(math.log10(distance[i,:,:,:].max())+1))
+
+# ------------------------------------------ assemble header ---------------------------------------
+
+ output = '%i\theader\n'%(len(new_header))
+ output += ''.join(new_header)
+
+# ------------------------------------- regenerate texture information ----------------------------------
+
+ for z in xrange(info['resolution'][2]):
+ for y in xrange(info['resolution'][1]):
+ output += {True:' ',False:'\n'}[options.twoD].join(map(lambda x: ('%%%ii'%formatwidth)%(round(x)), distance[i,:,y,z])) + '\n'
+
+
+# ------------------------------------------ output result ---------------------------------------
+
+ file['output'][i].write(output)
+
+ if file['name'] != 'STDIN':
+ file['output'][i].close()
+
+ if file['name'] != 'STDIN':
+ file['input'].close() # close input geom file
+
+