From 73c6bd767fde547876bdebf93bb2274dc6153b7e Mon Sep 17 00:00:00 2001 From: Tias Maiti Date: Tue, 26 May 2015 20:13:35 +0000 Subject: [PATCH] =?UTF-8?q?added=20options=20for=20selective=20seed=20pick?= =?UTF-8?q?ing=20based=20on=20Mitchell=E2=80=99s=20best=20candidate=20algo?= =?UTF-8?q?rithm=20for=20more=20uniformly=20distributed=20(spatially)=20se?= =?UTF-8?q?eds=20points?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- processing/pre/seeds_fromRandom.py | 90 ++++++++++++++++++++++-------- 1 file changed, 68 insertions(+), 22 deletions(-) diff --git a/processing/pre/seeds_fromRandom.py b/processing/pre/seeds_fromRandom.py index 4d11a6933..5cc3403f6 100755 --- a/processing/pre/seeds_fromRandom.py +++ b/processing/pre/seeds_fromRandom.py @@ -3,8 +3,10 @@ import os,sys,string,math,random import numpy as np -from optparse import OptionParser import damask +from optparse import OptionParser,OptionGroup +from scipy import spatial + scriptID = string.replace('$Id$','\n','\\n') scriptName = os.path.splitext(scriptID.split()[1])[0] @@ -33,6 +35,17 @@ parser.add_option('--sigma', dest='sigma', type='float', metavar='float', \ help='standard deviation of Gaussian Distribution for weights [%default]') parser.add_option('-m', '--microstructure', dest='microstructure', type='int', help='first microstructure index [%default]', metavar='int') +parser.add_option('-s','--selective', dest='selective', action='store_true', + help = 'selective picking of seed points from random seed points [%default]') +group = OptionGroup(parser, "Selective Seeding Options", + "More uniform distribution of seed points using Mitchells Best Candidate Algorithm" + ) +group.add_option('--distance', dest='bestDistance', type='float', metavar='float', \ + help='minimum distance to the next neighbor [%default]') +group.add_option('--numCandidates', dest='numCandidates', type='int', metavar='int', \ + help='maximum number of point to consider for initial random points generation [%default]') +parser.add_option_group(group) + parser.set_defaults(randomSeed = None) parser.set_defaults(grid = (16,16,16)) parser.set_defaults(N = 20) @@ -40,6 +53,10 @@ parser.set_defaults(weights=False) parser.set_defaults(mean = 0.0) parser.set_defaults(sigma = 1.0) parser.set_defaults(microstructure = 1) +parser.set_defaults(selective = False) +parser.set_defaults(bestDistance = 0.2) +parser.set_defaults(numCandidates = 10) + (options,filename) = parser.parse_args() @@ -47,6 +64,18 @@ options.grid = np.array(options.grid) labels = "1_coords\t2_coords\t3_coords\tphi1\tPhi\tphi2\tmicrostructure" +# ------------------------------------------ Functions Definitions --------------------------------- + +def kdtree_search(xyz, point) : + dist, index = spatial.cKDTree(xyz).query(np.array(point)) + return dist + +def generatePoint() : + return np.array([random.uniform(0,options.grid[0]/max(options.grid)), \ + random.uniform(0,options.grid[1]/max(options.grid)), \ + random.uniform(0,options.grid[2]/max(options.grid))]) + + # ------------------------------------------ setup file handle ------------------------------------- if filename == []: file = {'output':sys.stdout, 'croak':sys.stderr} @@ -69,28 +98,45 @@ grainEuler[0,:] *= 360.0 grainEuler[1,:] = np.arccos(2*grainEuler[1,:]-1)*180.0/math.pi # cos(Phi) is uniformly distributed grainEuler[2,:] *= 360.0 # phi_2 is uniformly distributed -seedpoints = -np.ones(options.N,dtype='int') # init grid positions of seed points - -if options.N * 1024 < gridSize: # heuristic limit for random search - i = 0 - while i < options.N: # until all (unique) points determined - p = np.random.randint(gridSize) # pick a location - if p not in seedpoints: # not yet taken? - seedpoints[i] = p # take it - i += 1 # advance stepper -else: - seedpoints = np.array(random.sample(range(gridSize),options.N)) # create random permutation of all grid positions and choose first N - -seeds = np.zeros((3,options.N),float) # init seed positions -seeds[0,:] = (np.mod(seedpoints ,options.grid[0])\ - +np.random.random())/options.grid[0] -seeds[1,:] = (np.mod(seedpoints// options.grid[0] ,options.grid[1])\ - +np.random.random())/options.grid[1] -seeds[2,:] = (np.mod(seedpoints//(options.grid[1]*options.grid[0]),options.grid[2])\ - +np.random.random())/options.grid[2] microstructure=np.arange(options.microstructure,options.microstructure+options.N).reshape(1,options.N) -table = np.transpose(np.concatenate((seeds,grainEuler,microstructure),axis = 0)) +if options.selective == False : + seedpoints = -np.ones(options.N,dtype='int') # init grid positions of seed points + + if options.N * 1024 < gridSize: # heuristic limit for random search + i = 0 + while i < options.N: # until all (unique) points determined + p = np.random.randint(gridSize) # pick a location + if p not in seedpoints: # not yet taken? + seedpoints[i] = p # take it + i += 1 # advance stepper + else: + seedpoints = np.array(random.sample(range(gridSize),options.N)) # create random permutation of all grid positions and choose first N + + seeds = np.zeros((3,options.N),float) # init seed positions + seeds[0,:] = (np.mod(seedpoints ,options.grid[0])\ + +np.random.random())/options.grid[0] + seeds[1,:] = (np.mod(seedpoints// options.grid[0] ,options.grid[1])\ + +np.random.random())/options.grid[1] + seeds[2,:] = (np.mod(seedpoints//(options.grid[1]*options.grid[0]),options.grid[2])\ + +np.random.random())/options.grid[2] + table = np.transpose(np.concatenate((seeds,grainEuler,microstructure),axis = 0)) +else : + samples = generatePoint().reshape(1,3) + + while samples.shape[0] < options.N : + bestDistance = options.bestDistance + for i in xrange(options.numCandidates) : + c = generatePoint() + d = kdtree_search(samples, c) + if (d > bestDistance) : + bestDistance = d + bestCandidate = c + if kdtree_search(samples,bestCandidate) != 0.0 : + samples = np.append(samples,bestCandidate.reshape(1,3),axis=0) + else : + continue + table = np.transpose(np.concatenate((samples.T,grainEuler,microstructure),axis = 0)) if options.weights : weight = np.random.normal(loc=options.mean, scale=options.sigma, size=options.N) @@ -98,7 +144,7 @@ if options.weights : table = np.append(table, weight.reshape(options.N,1), axis=1) labels += "\tweight" - +# -------------------------------------- Write Data -------------------------------------------------- header = ["5\theader", scriptID + " " + " ".join(sys.argv[1:]),