diff --git a/processing/pre/OIMlinear2linearODF.py b/processing/pre/OIMlinear2linearODF.py new file mode 100755 index 000000000..5277a8e1c --- /dev/null +++ b/processing/pre/OIMlinear2linearODF.py @@ -0,0 +1,67 @@ +#!/usr/bin/env python + +import os,sys,re + +try: + inName = sys.argv[1] + outName = os.path.splitext(inName)[0]+'.linearODF' + nPhi1,nPHI,nPhi2 = map(int,sys.argv[2:5]) +except: + print "\nusage:",sys.argv[0],"file nPhi1 nPHI nPhi2\n" + sys.exit(1) + +N = (nPhi1-1)*(nPHI-1)*(nPhi2-1) + +try: + inFile = open(inName,'r') + content = inFile.readlines() +except: + print 'unable to read:',inName + sys.exit(1) +try: + outFile = open(outName,'w') +except: + print 'unable to write:',outName + sys.exit(1) + +ODF = [[[[None] for k in range(nPhi2)] for j in range(nPHI)] for i in range(nPhi1)] +linear = [None]*N +line = 0 + +while (content[line].startswith('#')): # skip comments at start of file + line += 1 + +for iPhi1 in range(nPhi1): + for iPHI in range(nPHI): + for iPhi2 in range(nPhi2): + words = content[line].split() + ODF[iPhi1][iPHI][iPhi2] = float(words[3]) # extract intensity (in column 4) + line += 1 + +for iPhi1 in range(nPhi1-1): + for iPHI in range(nPHI-1): + for iPhi2 in range(nPhi2-1): + linear[iPhi1*(nPHI-1)*(nPhi2-1)+iPHI*(nPhi2-1)+iPhi2] = (\ + ODF[iPhi1 ][iPHI ][iPhi2 ] + \ + ODF[iPhi1 ][iPHI ][iPhi2+1] + \ + ODF[iPhi1 ][iPHI+1][iPhi2 ] + \ + ODF[iPhi1 ][iPHI+1][iPhi2+1] + \ + ODF[iPhi1+1][iPHI ][iPhi2 ] + \ + ODF[iPhi1+1][iPHI ][iPhi2+1] + \ + ODF[iPhi1+1][iPHI+1][iPhi2 ] + \ + ODF[iPhi1+1][iPHI+1][iPhi2+1] \ + ) / 8.0 + + +inFile.close() + +outFile.write('rangePhi1?\trangePHI?\trangePhi2?\n') +outFile.write('%i\t%i\t%i needs to be converted to angular steps\n'%(nPhi1-1,nPHI-1,nPhi2-1)) +outFile.write('cell-centered data\n') +outFile.write('\n') + +for i in range(N): + outFile.write('%g\n'%(linear[i])) + +outFile.close() + diff --git a/processing/pre/cropLinearODF.py b/processing/pre/cropLinearODF.py new file mode 100755 index 000000000..f7f695c1c --- /dev/null +++ b/processing/pre/cropLinearODF.py @@ -0,0 +1,105 @@ +#!/usr/bin/env python + +import os,sys,math,re + +# --- helper functions --- + +def binAsBins(bin,intervals): + """ explode compound bin into 3D bins list """ + bins = [0]*3 + bins[0] = (bin//(intervals[1] * intervals[2])) % intervals[0] + bins[1] = (bin//intervals[2]) % intervals[1] + bins[2] = bin % intervals[2] + return bins + +def binsAsBin(bins,intervals): + """ implode 3D bins into compound bin """ + return (bins[0]*intervals[1] + bins[1])*intervals[2] + bins[2] + +def EulersAsBins(Eulers,intervals,deltas,center): + """ return list of Eulers translated into 3D bins list """ + return [\ + int((euler+(0.5-center)*delta)//delta)%interval \ + for euler,delta,interval in zip(Eulers,deltas,intervals) \ + ] + +def binAsEulers(bin,intervals,deltas,center): + """ compound bin number translated into list of Eulers """ + Eulers = [0.0]*3 + Eulers[2] = (bin%intervals[2] + center)*deltas[2] + Eulers[1] = (bin//intervals[2]%intervals[1] + center)*deltas[1] + Eulers[0] = (bin//(intervals[2]*intervals[1]) + center)*deltas[0] + return Eulers + + + +# check usage +try: + inName = sys.argv[1] + outLimits = sys.argv[2:5] +except: + print "usage:",sys.argv[0],"nameLinearODF limitPhi1 limitPHI limitPhi2" + sys.exit(1); + + +#open binned ODF +try: + inFile = open(inName,'r') +except: + print 'unable to open binnedODF:', inName; + sys.exit(1); + +# process header info +ODF = {} +inLimits = [math.radians(int(float(limit))) for limit in inFile.readline().split()] +outLimits = [math.radians(int(float(limit))) for limit in outLimits] +deltas = [math.radians(float(delta)) for delta in inFile.readline().split()] +inIntervals = [int(limit/delta) for limit,delta in zip(inLimits,deltas)] +outIntervals = [int(limit/delta) for limit,delta in zip(outLimits,deltas)] +inBins = inIntervals[0]*inIntervals[1]*inIntervals[2] + +print 'Limit:', [math.degrees(limit) for limit in inLimits] +print 'Crop:', [math.degrees(limit) for limit in outLimits] +print 'Delta:', [math.degrees(delta) for delta in deltas] +print 'Interval:', inIntervals +print 'Interval:', outIntervals + +centering = inFile.readline() +if re.search('cell',centering.lower()): + ODF['center'] = 0.5 + print 'cell-centered data (offset %g)'%ODF['center'] +else: + ODF['center'] = 0.0 + print 'vertex-centered data (offset %g)'%ODF['center'] + +inFile.readline() # skip blank delimiter + +# read linear binned data +inODF = map(float,inFile.readlines()) +inFile.close() + +if len(inODF) != inBins: + print 'expecting', inBins, 'values but got', len(inODF) + sys.exit(1) + +try: + outName = os.path.splitext(inName)[0]+'_%ix%ix%i'%(outIntervals[0],outIntervals[1],outIntervals[2])+'.linearODF' + outFile = open(outName,'w') +except: + print 'unable to write:',outName + sys.exit(1) + +outFile.write('%g\t%g\t%g\n'%(\ + math.degrees(outIntervals[0]*deltas[0]),\ + math.degrees(outIntervals[1]*deltas[1]),\ + math.degrees(outIntervals[2]*deltas[2]) )) +outFile.write('%i\t%i\t%i\n'%(math.degrees(deltas[0]),math.degrees(deltas[1]),math.degrees(deltas[2]))) +outFile.write('%s-centered data\n'%{True:'vertex',False:'cell'}[ODF['center']==0.0]) +outFile.write('\n') + +for phi1 in range(outIntervals[0]): + for Phi in range(outIntervals[1]): + for phi2 in range(outIntervals[2]): + outFile.write('%g\n'%(inODF[((phi1*inIntervals[1])+Phi)*inIntervals[2]+phi2])) + +outFile.close() diff --git a/processing/pre/hybridIA_linODFsampling.py b/processing/pre/hybridIA_linODFsampling.py new file mode 100755 index 000000000..b53989284 --- /dev/null +++ b/processing/pre/hybridIA_linODFsampling.py @@ -0,0 +1,331 @@ +#!/usr/bin/env python + +import os,sys,math,re,random +random.seed() + + +# --- helper functions --- + +def binAsBins(bin,intervals): + """ explode compound bin into 3D bins list """ + bins = [0]*3 + bins[0] = (bin//(intervals[1] * intervals[2])) % intervals[0] + bins[1] = (bin//intervals[2]) % intervals[1] + bins[2] = bin % intervals[2] + return bins + +def binsAsBin(bins,intervals): + """ implode 3D bins into compound bin """ + return (bins[0]*intervals[1] + bins[1])*intervals[2] + bins[2] + +def EulersAsBins(Eulers,intervals,deltas,center): + """ return list of Eulers translated into 3D bins list """ + return [\ + int((euler+(0.5-center)*delta)//delta)%interval \ + for euler,delta,interval in zip(Eulers,deltas,intervals) \ + ] + +def binAsEulers(bin,intervals,deltas,center): + """ compound bin number translated into list of Eulers """ + Eulers = [0.0]*3 + Eulers[2] = (bin%intervals[2] + center)*deltas[2] + Eulers[1] = (bin//intervals[2]%intervals[1] + center)*deltas[1] + Eulers[0] = (bin//(intervals[2]*intervals[1]) + center)*deltas[0] + return Eulers + +def directInvRepetitions(probability,scale): + """ calculate number of samples drawn by direct inversion """ + nDirectInv = 0 + for bin in range(len(probability)): # loop over bins + nDirectInv += int(round(probability[bin]*scale)) # calc repetition + return nDirectInv + + +# --- sampling methods --- + +# ----- efficient algorithm --------- + +def directInversion (ODF,nSamples): + """ ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians) """ + + nBins = ODF['intervals'][0]*ODF['intervals'][1]*ODF['intervals'][2] + deltas = [limit/intervals for limit,intervals in zip(ODF['limits'],ODF['intervals'])] + + # calculate repetitions of each orientation + if re.search(r'hybrid',sys.argv[0],re.IGNORECASE): # my script's name contains "hybrid" + nOptSamples = max(ODF['nNonZero'],nSamples) # random subsampling if too little samples requested + else: # blunt integer approximation + nOptSamples = nSamples + + nInvSamples = 0 + repetition = [None]*nBins + probabilityScale = nOptSamples # guess + + scaleLower = 0.0 + nInvSamplesLower = 0 + scaleUpper = float(nOptSamples) + incFactor = 1.0 + nIter = 0 + nInvSamplesUpper = directInvRepetitions(ODF['dV_V'],scaleUpper) + while (\ + (scaleUpper-scaleLower > scaleUpper*1e-15 or nInvSamplesUpper < nOptSamples) and \ + nInvSamplesUpper != nOptSamples \ + ): # closer match required? + if nInvSamplesUpper < nOptSamples: + scaleLower,scaleUpper = scaleUpper,scaleUpper+incFactor*(scaleUpper-scaleLower)/2.0 + incFactor *= 2.0 + nInvSamplesLower,nInvSamplesUpper = nInvSamplesUpper,directInvRepetitions(ODF['dV_V'],scaleUpper) + else: + scaleUpper = (scaleLower+scaleUpper)/2.0 + incFactor = 1.0 + nInvSamplesUpper = directInvRepetitions(ODF['dV_V'],scaleUpper) + nIter += 1 + print '%i:(%12.11f,%12.11f) %i <= %i <= %i'%(nIter,scaleLower,scaleUpper,nInvSamplesLower,nOptSamples,nInvSamplesUpper) + nInvSamples = nInvSamplesUpper + scale = scaleUpper + print 'created set of',nInvSamples,'samples (',float(nInvSamples)/nOptSamples-1.0,') with scaling',scale,'delivering',nSamples + repetition = [None]*nBins # preallocate and clear + + for bin in range(nBins): # loop over bins + repetition[bin] = int(round(ODF['dV_V'][bin]*scale)) # calc repetition + + # build set + set = [None]*nInvSamples + i = 0 + for bin in range(nBins): + set[i:i+repetition[bin]] = [bin]*repetition[bin] # fill set with bin, i.e. orientation + i += repetition[bin] # advance set counter + + orientations = [None]*nSamples + reconstructedODF = [0.0]*nBins + unitInc = 1.0/nSamples + for j in range(nSamples): + if (j == nInvSamples-1): ex = j + else: ex = int(round(random.uniform(j+0.5,nInvSamples-0.5))) + bin = set[ex] + bins = binAsBins(bin,ODF['intervals']) + Eulers = binAsEulers(bin,ODF['intervals'],deltas,ODF['center']) + orientations[j] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) ) + reconstructedODF[bin] += unitInc + set[ex] = set[j] # exchange orientations + + return orientations, reconstructedODF + + +# ----- trial and error algorithms --------- + +def MonteCarloEulers (ODF,nSamples): + """ ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians) """ + + countMC = 0 + maxdV_V = max(ODF['dV_V']) + nBins = ODF['intervals'][0]*ODF['intervals'][1]*ODF['intervals'][2] + deltas = [limit/intervals for limit,intervals in zip(ODF['limits'],ODF['intervals'])] + orientations = [None]*nSamples + reconstructedODF = [0.0]*nBins + unitInc = 1.0/nSamples + + for j in range(nSamples): + MC = maxdV_V*2.0 + bin = 0 + while MC > ODF['dV_V'][bin]: + countMC += 1 + MC = maxdV_V*random.random() + Eulers = [limit*random.random() for limit in ODF['limits']] + bins = EulersAsBins(Eulers,ODF['intervals'],deltas,ODF['center']) + bin = binsAsBin(bins,ODF['intervals']) + orientations[j] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) ) + reconstructedODF[bin] += unitInc + + return orientations, reconstructedODF, countMC + + +def MonteCarloBins (ODF,nSamples): + """ ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians) """ + + countMC = 0 + maxdV_V = max(ODF['dV_V']) + nBins = ODF['intervals'][0]*ODF['intervals'][1]*ODF['intervals'][2] + deltas = [limit/intervals for limit,intervals in zip(ODF['limits'],ODF['intervals'])] + orientations = [None]*nSamples + reconstructedODF = [0.0]*nBins + unitInc = 1.0/nSamples + + for j in range(nSamples): + MC = maxdV_V*2.0 + bin = 0 + while MC > ODF['dV_V'][bin]: + countMC += 1 + MC = maxdV_V*random.random() + bin = int(nBins * random.random()) + Eulers = binAsEulers(bin,ODF['intervals'],deltas,ODF['center']) + orientations[j] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) ) + reconstructedODF[bin] += unitInc + + return orientations, reconstructedODF + + +def TothVanHoutteSTAT (ODF,nSamples): + """ ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians) """ + + nBins = ODF['intervals'][0]*ODF['intervals'][1]*ODF['intervals'][2] + deltas = [limit/intervals for limit,intervals in zip(ODF['limits'],ODF['intervals'])] + orientations = [None]*nSamples + reconstructedODF = [0.0]*nBins + unitInc = 1.0/nSamples + + selectors = [random.random() for i in range(nSamples)] + selectors.sort() + indexSelector = 0 + + cumdV_V = 0.0 + countSamples = 0 + + for bin in range(nBins) : + cumdV_V += ODF['dV_V'][bin] + while indexSelector < nSamples and selectors[indexSelector] < cumdV_V: + Eulers = binAsEulers(bin,ODF['intervals'],deltas,ODF['center']) + orientations[countSamples] = '%g\t%g\t%g' %( math.degrees(Eulers[0]),math.degrees(Eulers[1]),math.degrees(Eulers[2]) ) + reconstructedODF[bin] += unitInc + countSamples += 1 + indexSelector += 1 + + print 'created set of',countSamples,'when asked to deliver',nSamples + + return orientations, reconstructedODF + + + + +# check usage +try: + nameBinnedODF = sys.argv[1] + nSamples = int(float(sys.argv[2])) +except: + print "\nusage:",os.path.basename(sys.argv[0]),"nameLinearODF nSamples [nameSampledODF]\n" + sys.exit(1); + +methods = ['IA','STAT'] + +#open binned ODF +try: + fileBinnedODF = open(nameBinnedODF,'r') +except: + print 'unable to open binnedODF:', nameBinnedODF; + sys.exit(1); + +# process header info +ODF = {} +ODF['limits'] = [math.radians(float(limit)) for limit in fileBinnedODF.readline().split()] +ODF['deltas'] = [math.radians(float(delta)) for delta in fileBinnedODF.readline().split()] +ODF['intervals'] = [int(round(limit/delta)) for limit,delta in zip(ODF['limits'],ODF['deltas'])] +nBins = ODF['intervals'][0]*ODF['intervals'][1]*ODF['intervals'][2] + +print 'Limit:', [math.degrees(limit) for limit in ODF['limits']] +print 'Delta:', [math.degrees(delta) for delta in ODF['deltas']] +print 'Interval:', ODF['intervals'] + +centering = fileBinnedODF.readline() +if re.search('cell',centering.lower()): + ODF['center'] = 0.5 + print 'cell-centered data (offset %g)'%ODF['center'] +else: + ODF['center'] = 0.0 + print 'vertex-centered data (offset %g)'%ODF['center'] + +fileBinnedODF.readline() # skip blank delimiter + +# read linear binned data +linesBinnedODF = fileBinnedODF.readlines() +fileBinnedODF.close() + +if len(linesBinnedODF) != nBins: + print 'expecting', nBins, 'values but got', len(linesBinnedODF) + sys.exit(1) + +# build binnedODF array +print 'reading',nBins,'values' +sumdV_V = 0.0 +ODF['dV_V'] = [None]*nBins +ODF['nNonZero'] = 0 +dg = ODF['deltas'][0]*2*math.sin(ODF['deltas'][1]/2.0)*ODF['deltas'][2] +for bin in range(nBins): + ODF['dV_V'][bin] = \ + max(0.0,float(linesBinnedODF[bin])) * dg * \ + math.sin(((bin//ODF['intervals'][2])%ODF['intervals'][1]+ODF['center'])*ODF['deltas'][1]) + if ODF['dV_V'][bin] > 0.0: + sumdV_V += ODF['dV_V'][bin] + ODF['nNonZero'] += 1 + +for bin in range(nBins): ODF['dV_V'][bin] /= sumdV_V # normalize dV/V + +print 'non-zero fraction:', float(ODF['nNonZero'])/nBins,'(%i/%i)'%(ODF['nNonZero'],nBins) +print 'Volume integral of ODF:', sumdV_V +print 'Reference Integral:', ODF['limits'][0]*ODF['limits'][2]*(1-math.cos(ODF['limits'][1])) + +# call methods +Functions = {'IA': 'directInversion', 'STAT': 'TothVanHoutteSTAT', 'MC': 'MonteCarloBins'} +Orientations = {} +ReconstructedODF = {} +for method in methods: + Orientations[method], ReconstructedODF[method] = (globals()[Functions[method]])(ODF,nSamples) + +# calculate accuracy of sample +squaredDiff = {} +squaredRelDiff = {} +mutualProd = {} +indivSum = {} +indivSquaredSum = {} +for method in ['orig']+methods: + squaredDiff[method] = 0.0 + squaredRelDiff[method] = 0.0 + mutualProd[method] = 0.0 + indivSum[method] = 0.0 + indivSquaredSum[method] = 0.0 + +for bin in range(nBins): + for method in methods: + squaredDiff[method] += (ODF['dV_V'][bin] - ReconstructedODF[method][bin])**2 + if ODF['dV_V'][bin] > 0.0: + squaredRelDiff[method] += (ODF['dV_V'][bin] - ReconstructedODF[method][bin])**2/ODF['dV_V'][bin]**2 + mutualProd[method] += ODF['dV_V'][bin]*ReconstructedODF[method][bin] + indivSum[method] += ReconstructedODF[method][bin] + indivSquaredSum[method] += ReconstructedODF[method][bin]**2 + indivSum['orig'] += ODF['dV_V'][bin] + indivSquaredSum['orig'] += ODF['dV_V'][bin]**2 + +print 'sqrt(N*)RMSD of ODFs:\t', [math.sqrt(nSamples*squaredDiff[method]) for method in methods] +print 'RMSrD of ODFs:\t', [math.sqrt(squaredRelDiff[method]) for method in methods] +print 'rMSD of ODFs:\t', [squaredDiff[method]/indivSquaredSum['orig'] for method in methods] +#print 'correlation slope:\t', [(nBins*mutualProd[method]-indivSum['orig']*indivSum[method])/(nBins*indivSquaredSum['orig']-indivSum['orig']**2) for method in ('IA','STAT','MC')] +#print 'correlation confidence:\t', [(mutualProd[method]-indivSum['orig']*indivSum[method]/nBins)/\ +# (nBins*math.sqrt((indivSquaredSum['orig']/nBins-(indivSum['orig']/nBins)**2)*(indivSquaredSum[method]/nBins-(indivSum[method]/nBins)**2))) for method in ('IA','STAT','MC')] +print 'nNonZero correlation slope:\t', [(ODF['nNonZero']*mutualProd[method]-indivSum['orig']*indivSum[method])/(ODF['nNonZero']*indivSquaredSum['orig']-indivSum['orig']**2) for method in methods] +print 'nNonZero correlation confidence:\t', [(mutualProd[method]-indivSum['orig']*indivSum[method]/ODF['nNonZero'])/\ + (ODF['nNonZero']*math.sqrt((indivSquaredSum['orig']/ODF['nNonZero']-(indivSum['orig']/ODF['nNonZero'])**2)*(indivSquaredSum[method]/ODF['nNonZero']-(indivSum[method]/ODF['nNonZero'])**2))) for method in methods] + +# write result +try: + nameSampledODF = sys.argv[3] +except: + sys.exit(0) # that's it folks + +for method in methods: + if method == 'IA' and nSamples < ODF['nNonZero']: + strOpt = '(%i)'%ODF['nNonZero'] + else: + strOpt = '' + try: + fileSampledODF = open(nameSampledODF+'.'+method+'sampled_'+str(nSamples)+strOpt, 'w') + fileSampledODF.write('%i\n'%nSamples) + fileSampledODF.write('\n'.join(Orientations[method])+'\n') + fileSampledODF.close() + except: + print 'unable to write sampledODF:', nameSampledODF+'.'+method+'sampled_'+str(nSamples)+strOpt + + try: + fileRegressionODF = open(nameSampledODF+'.'+method+'regression_'+str(nSamples)+strOpt, 'w') + fileRegressionODF.write('\n'.join([a+'\t'+b for (a,b) in zip(map(str,ReconstructedODF[method]),map(str,ODF['dV_V']))])+'\n') + fileRegressionODF.close() + except: + print 'unable to write RegressionODF:', nameSampledODF+'.'+method+'regression_'+str(nSamples)+strOpt diff --git a/processing/pre/texture2ang.py b/processing/pre/texture2ang.py new file mode 100755 index 000000000..8697ece5e --- /dev/null +++ b/processing/pre/texture2ang.py @@ -0,0 +1,121 @@ +#!/usr/bin/env python + +import os,sys,math,re +from optparse import OptionParser + +def integerFactorization(i): + + j = int(math.floor(math.sqrt(float(i)))) + while (j>1 and int(i)%j != 0): + j -= 1 + return j + +def positiveRadians(angle): + + angle = math.radians(float(angle)) + while angle < 0.0: + angle += 2.0*math.pi + + return angle + + +def getHeader(sizeX,sizeY,step): + + return [ \ + '# TEM_PIXperUM 1.000000', \ + '# x-star 0.509548', \ + '# y-star 0.795272', \ + '# z-star 0.611799', \ + '# WorkingDistance 18.000000', \ + '#', \ + '# Phase 1', \ + '# MaterialName Al', \ + '# Formula Fe', \ + '# Info', \ + '# Symmetry 43', \ + '# LatticeConstants 2.870 2.870 2.870 90.000 90.000 90.000', \ + '# NumberFamilies 4', \ + '# hklFamilies 1 1 0 1 0.000000 1', \ + '# hklFamilies 2 0 0 1 0.000000 1', \ + '# hklFamilies 2 1 1 1 0.000000 1', \ + '# hklFamilies 3 1 0 1 0.000000 1', \ + '# Categories 0 0 0 0 0 ', \ + '#', \ + '# GRID: SquareGrid', \ + '# XSTEP: ' + str(step), \ + '# YSTEP: ' + str(step), \ + '# NCOLS_ODD: ' + str(sizeX), \ + '# NCOLS_EVEN: ' + str(sizeX), \ + '# NROWS: ' + str(sizeY), \ + '#', \ + '# OPERATOR: ODFsammpling', \ + '#', \ + '# SAMPLEID: ', \ + '#', \ + '# SCANID: ', \ + '#', \ + ] + + +parser = OptionParser(usage='%prog [options] datafile(s)') +parser.add_option("-c", "--column", type="int",\ + dest="column",\ + help="starting column of Euler triplet") +parser.add_option("-s", "--skip", type="int",\ + dest="skip",\ + help="skip this many lines of heading info [%default]") + +parser.set_defaults (column = 1) +parser.set_defaults (skip = 0) + +(options, files) = parser.parse_args() +options.column -= 1 + +if files == []: + parser.error('no input file specified...') + sys.exit(1) +if options.column < 0: + parser.error('column needs to be 1,2,...') + sys.exit(1) + +while len(files) > 0: + textureFilename = files.pop() + baseName = os.path.splitext(textureFilename)[0] + + # open texture file and read content + textureFile = open(textureFilename) + content = textureFile.readlines() + textureFile.close() + + m = re.match('(\d+)\s+head',content[0],re.I) + if m != None and options.skip == 0: + options.skip = int(m.group(1))+1 + + # extract orientation angles + angles = [map(positiveRadians,line.split()[options.column:options.column+3]) for line in content[options.skip:]] + + nPoints = len(angles) + sizeY = integerFactorization(nPoints) + sizeX = nPoints / sizeY + print '%s: %i * %i = %i (== %i)'%(baseName,sizeX,sizeY,sizeX*sizeY,nPoints) + # write ang file + + try: + + # write header + angFile = open(baseName + '.ang','w') + for line in getHeader(sizeX,sizeY,1.0): + angFile.write(line + '\n') + + # write data + counter = 0 + for point in angles: + angFile.write(''.join(['%10.5f'%angle for angle in point])+ + ''.join(['%10.5f'%coord for coord in [counter%sizeX,counter//sizeX]])+ + ' 100.0 1.0 0 1 1.0\n') + counter += 1 + + angFile.close() + + except: + print 'unable to write',baseName