Skip to content

Commit b6cae6c

Browse files
author
Joel Bernier
committed
added overlap table script
1 parent ee844bf commit b6cae6c

File tree

2 files changed

+228
-8
lines changed

2 files changed

+228
-8
lines changed

‎hexrd/fitgrains.py‎

Lines changed: 1 addition & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -344,7 +344,7 @@ def fit_grains(self, grain_id, grain_params, refit_tol=None):
344344
#print "no overlap table found"
345345
pass
346346

347-
# completeness from pullspots only; incl saturated
347+
# completeness from pullspots only; incl saturated and overlaps
348348
completeness = sum(valid_refl_ids)/float(len(valid_refl_ids))
349349

350350
# extract data from grain table
@@ -447,13 +447,6 @@ def get_residuals(self, grain_params):
447447

448448
def loop(self):
449449
id, grain_params = self._jobs.get(False)
450-
451-
# slated for removal ## skips the first loop if have_estimate is True
452-
# slated for removal #if self._p['skip_on_estimate']:
453-
# slated for removal # have_estimate = not np.all(grain_params[-9:] == [0,0,0,1,1,1,0,0,0])
454-
# slated for removal # iterations = (have_estimate, len(self._p['eta_tol']))
455-
# slated for removal #else:
456-
# slated for removal # iterations = (0, len(self._p['eta_tol']))
457450
iterations = (0, len(self._p['eta_tol']))
458451
for iteration in range(*iterations):
459452
# pull spots if asked to, otherwise just fit

‎scripts/makeOverlapTable.py‎

Lines changed: 227 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,227 @@
1+
# -*- coding: utf-8 -*-
2+
"""
3+
Spyder Editor
4+
5+
This is a temporary script file.
6+
"""
7+
8+
import argparse, os, sys
9+
10+
import cPickle
11+
12+
import numpy as np
13+
14+
from mpl_toolkits.mplot3d import Axes3D
15+
import matplotlib.pyplot as plt
16+
from matplotlib.colors import BoundaryNorm
17+
18+
from hexrd import config
19+
from hexrd.xrd import transforms_CAPI as xfcapi
20+
from hexrd.coreutil import get_instrument_parameters
21+
22+
from sklearn.cluster import dbscan
23+
24+
from scipy import cluster
25+
26+
def adist(ang0, ang1):
27+
resd = xfcapi.angularDifference(ang0 ,ang1)
28+
return np.sqrt(sum(resd**2))
29+
30+
def build_overlap_table(cfg, tol_mult=0.5):
31+
32+
icfg = get_instrument_parameters(cfg)
33+
34+
gt = np.loadtxt(
35+
os.path.join(cfg.analysis_dir, 'grains.out')
36+
)
37+
38+
ngrains = len(gt)
39+
40+
mat_list = cPickle.load(open(cfg.material.definitions, 'r'))
41+
mat_names = [mat_list[i].name for i in range(len(mat_list))]
42+
mat_dict = dict(zip(mat_names, mat_list))
43+
44+
matl = mat_dict[cfg.material.active]
45+
46+
pd = matl.planeData
47+
pd.exclusions = np.zeros(len(pd.exclusions), dtype=bool)
48+
pd.tThMax = np.radians(cfg.fit_grains.tth_max)
49+
pd.tThWidth = np.radians(cfg.fit_grains.tolerance.tth[-1])
50+
51+
# for clustering...
52+
eps = tol_mult*np.radians(
53+
min(
54+
min(cfg.fit_grains.tolerance.eta),
55+
2*min(cfg.fit_grains.tolerance.omega)
56+
)
57+
)
58+
59+
# merged two-theta indices
60+
tth_ranges_merged = pd.getMergedRanges()[0]
61+
pids = []
62+
for hklids in tth_ranges_merged:
63+
pids.append(
64+
[pd.hklDataList[hklids[i]]['hklID'] for i in range(len(hklids))]
65+
)
66+
67+
# Make table of unit diffraction vectors
68+
st = []
69+
for i in range(ngrains):
70+
this_st = np.loadtxt(
71+
os.path.join(cfg.analysis_dir, 'spots_%05d.out' %i)
72+
)
73+
#... do all predicted?
74+
valid_spt = this_st[:, 0] >= 0
75+
#valid_spt = np.ones(len(this_st), dtype=bool)
76+
77+
angs = this_st[valid_spt, 7:10]
78+
79+
dvec = xfcapi.anglesToDVec(
80+
angs,
81+
chi=icfg['oscillation_stage']['chi']
82+
)
83+
84+
# [ grainID, reflID, hklID, D_s[0], D_s[1], D_s[2], tth, eta, ome ]
85+
st.append(
86+
np.hstack([
87+
i*np.ones((sum(valid_spt), 1)),
88+
this_st[valid_spt, :2],
89+
dvec,
90+
angs,
91+
])
92+
)
93+
94+
# make overlap table
95+
# [[range_0], [range_1], ..., [range_n]]
96+
# range_0 = [grainIDs, reflIDs, hklIDs] that are within tol
97+
overlap_table = []
98+
ii = 0
99+
for pid in pids:
100+
tmp = []; a = []; b = []; c = []
101+
for j in range(len(pid)):
102+
a.append(
103+
np.vstack(
104+
[st[i][st[i][:, 2] == pid[j], 3:6] for i in range(len(st))]
105+
)
106+
)
107+
b.append(
108+
np.vstack(
109+
[st[i][st[i][:, 2] == pid[j], 0:3] for i in range(len(st))]
110+
)
111+
)
112+
c.append(
113+
np.vstack(
114+
[st[i][st[i][:, 2] == pid[j], 6:9] for i in range(len(st))]
115+
)
116+
)
117+
pass
118+
a = np.vstack(a)
119+
b = np.vstack(b)
120+
c = np.vstack(c)
121+
if len(a) > 0:
122+
# run dbscan
123+
core_samples, labels = dbscan(
124+
a,
125+
eps=eps,
126+
min_samples=2,
127+
metric='minkowski', p=2,
128+
)
129+
130+
cl = np.array(labels, dtype=int) # convert to array
131+
noise_points = cl == -1 # index for marking noise
132+
cl += 1 # move index to 1-based instead of 0
133+
cl[noise_points] = -1 # re-mark noise as -1
134+
135+
# extract number of clusters
136+
if np.any(cl == -1):
137+
nblobs = len(np.unique(cl)) - 1
138+
else:
139+
nblobs = len(np.unique(cl))
140+
141+
for i in range(1, nblobs+1):
142+
# put in check on omega here
143+
these_angs = c[np.where(cl == i)[0], :]
144+
local_cl = cluster.hierarchy.fclusterdata(
145+
these_angs[:, 1:],
146+
eps,
147+
criterion='distance',
148+
metric=adist
149+
)
150+
local_nblobs = len(np.unique(local_cl))
151+
if local_nblobs < len(these_angs):
152+
for j in range(1, local_nblobs + 1):
153+
npts = sum(local_cl == j)
154+
if npts >= 2:
155+
cl_idx = np.where(local_cl == j)[0]
156+
#import pdb; pdb.set_trace()
157+
tmp.append(
158+
b[np.where(cl == i)[0][cl_idx], :]
159+
)
160+
print "processing ring set %d" %ii
161+
ii += 1
162+
overlap_table.append(tmp)
163+
return overlap_table
164+
165+
def build_discrete_cmap(ngrains):
166+
167+
# define the colormap
168+
cmap = plt.cm.jet
169+
170+
# extract all colors from the .jet map
171+
cmaplist = [cmap(i) for i in range(cmap.N)]
172+
173+
# create the new map
174+
cmap = cmap.from_list('Custom cmap', cmaplist, cmap.N)
175+
176+
# define the bins and normalize
177+
bounds = np.linspace(0, ngrains, ngrains+1)
178+
norm = BoundaryNorm(bounds, cmap.N)
179+
180+
return cmap, norm
181+
182+
#%%
183+
if __name__ == '__main__':
184+
parser = argparse.ArgumentParser(description='Make overlap table from cfg file')
185+
parser.add_argument(
186+
'cfg', metavar='cfg_filename',
187+
type=str, help='a YAML config filename')
188+
parser.add_argument(
189+
'-m', '--multiplier',
190+
help='multiplier on angular tolerance',
191+
type=float, default=0.5)
192+
193+
args = vars(parser.parse_args(sys.argv[1:]))
194+
195+
cfg = config.open(args['cfg'])[0]
196+
print "loaded config file %s" %args['cfg']
197+
overlap_table = build_overlap_table(cfg)
198+
np.savez(os.path.join(cfg.analysis_dir, 'overlap_table.npz'), *overlap_table)
199+
#%%
200+
#fig = plt.figure()
201+
#ax = fig.add_subplot(111, projection='3d')
202+
#
203+
#etas = np.radians(np.linspace(0, 359, num=360))
204+
#cx = np.cos(etas)
205+
#cy = np.sin(etas)
206+
#cz = np.zeros_like(etas)
207+
#
208+
#ax.plot(cx, cy, cz, c='b')
209+
#ax.plot(cx, cz, cy, c='g')
210+
#ax.plot(cz, cx, cy, c='r')
211+
#ax.scatter3D(a[:, 0], a[:, 1], a[:, 2], c=b[:, 0], cmap=cmap, norm=norm, marker='o', s=20)
212+
#
213+
#ax.set_xlabel(r'$\mathbf{\mathrm{X}}_s$')
214+
#ax.set_ylabel(r'$\mathbf{\mathrm{Y}}_s$')
215+
#ax.set_zlabel(r'$\mathbf{\mathrm{Z}}_s$')
216+
#
217+
#ax.elev = 124
218+
#ax.azim = -90
219+
#
220+
#ax.axis('equal')
221+
#
222+
##fname = "overlaps_%03d.png"
223+
##for i in range(360):
224+
## ax.azim += i
225+
## fig.savefig(
226+
## fname %i, dpi=200, facecolor='w', edgecolor='w',
227+
## orientation='landcape')

0 commit comments

Comments
 (0)