Skip to content

Commit

Permalink
Merge pull request #1244 from psavery/migrate-to-hexrd-powder-calibrator
Browse files Browse the repository at this point in the history 
Migrate PowderCalibrator to the one in hexrd
  • Loading branch information
@joelvbernier
joelvbernier committed Jun 15, 2022
2 parents eb45995 + 5307626 commit b6289b3
Show file tree
Hide file tree
Showing 2 changed files with 81 additions and 217 deletions.
17 changes: 8 additions & 9 deletions hexrd/ui/calibration/calibration_runner.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -373,13 +373,15 @@ def run_calibration(self):
flags = HexrdConfig().get_statuses_instrument_format()
instr.calibration_flags = flags

instr_calibrator = run_calibration(picks, instr, materials)
img_dict = HexrdConfig().masked_images_dict

instr_calibrator = run_calibration(picks, instr, img_dict, materials)
self.write_instrument_to_hexrd_config(instr)

# Update the lattice parameters and overlays
overlays = self.active_overlays
for overlay, calibrator in zip(overlays, instr_calibrator.calibrators):
if calibrator.calibrator_type == 'powder':
if isinstance(calibrator, PowderCalibrator):
if calibrator.params.size == 0:
continue

Expand All @@ -389,7 +391,7 @@ def run_calibration(self):
HexrdConfig().flag_overlay_updates_for_material(mat_name)
if mat is HexrdConfig().active_material:
HexrdConfig().active_material_modified.emit()
elif calibrator.calibrator_type == 'laue':
else:
overlay.crystal_params = calibrator.params

# In case any overlays changed
Expand Down Expand Up @@ -708,8 +710,7 @@ def auto_pick_powder_points(self):
options = HexrdConfig().config['calibration']['powder']
self.instr = create_hedm_instrument()

# Assume there is only one image in each image series for now...
img_dict = {k: x[0] for k, x in HexrdConfig().imageseries_dict.items()}
img_dict = HexrdConfig().masked_images_dict

statuses = HexrdConfig().get_statuses_instrument_format()
self.instr.calibration_flags = statuses
Expand Down Expand Up @@ -787,14 +788,12 @@ def auto_pick_laue_spots(self):
self.async_runner.run(self.run_auto_laue_pick)

def run_auto_laue_pick(self):
# Assume there is only one image in each image series for now...
imsd = HexrdConfig().imageseries_dict
raw_img_dict = {k: x[0] for k, x in imsd.items()}
img_dict = HexrdConfig().masked_images_dict

# These are the options the user chose earlier...
options = HexrdConfig().config['calibration']['laue_auto_picker']
kwargs = {
'raw_img_dict': raw_img_dict,
'raw_img_dict': img_dict,
**options
}
return self.laue_auto_picker._autopick_points(**kwargs)
Expand Down
281 changes: 73 additions & 208 deletions hexrd/ui/calibration/pick_based_calibration.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,9 @@
gaussian_2d, gaussian_2d_int, sxcal_obj_func,
__reflInfo_dtype as reflInfo_dtype
)
from hexrd.ui.utils.conversions import angles_to_cart
from hexrd.constants import fwhm_to_sigma
from hexrd.fitting.calibration import PowderCalibrator
from hexrd.transforms import xfcapi
from hexrd import xrdutil

Expand Down Expand Up @@ -479,202 +481,10 @@ def model(self, reduced_params, data_dict):
)


class PowderCalibrator(object):
_CALIBRATOR_TYPE = 'powder'

def __init__(self, instr, plane_data, flags):
self._instr = instr
self._plane_data = plane_data
self._plane_data.wavelength = self._instr.beam_energy # force
self._params = np.asarray(self._plane_data.lparms, dtype=float)
self._full_params = np.hstack(
[self._instr.calibration_parameters, self._params]
)
assert len(flags) == len(self._full_params), \
"flags must have %d elements" % len(self._full_params)
self._flags = flags

@property
def calibrator_type(self):
return self._CALIBRATOR_TYPE

@property
def instr(self):
return self._instr

@property
def plane_data(self):
self._plane_data.wavelength = self._instr.beam_energy
return self._plane_data

@property
def params(self):
return self._params

@params.setter
def params(self, x):
x = np.atleast_1d(x)
if len(x) != len(self.plane_data.lparms):
raise RuntimeError("params must have %d elements"
% len(self.plane_data.lparms))
self._params = x
self._plane_data.lparms = x

@property
def full_params(self):
return self._full_params

@property
def npi(self):
return len(self._instr.calibration_parameters)

@property
def npe(self):
return len(self._params)

@property
def flags(self):
return self._flags

@flags.setter
def flags(self, x):
x = np.atleast_1d(x)
nparams_instr = len(self.instr.calibration_parameters)
nparams_extra = len(self.params)
nparams = nparams_instr + nparams_extra
if len(x) != nparams:
raise RuntimeError("flags must have %d elements" % nparams)
self._flags = np.asarrasy(x, dtype=bool)
self._instr.calibration_flags = self._flags[:nparams_instr]

def _evaluate(self, reduced_params, data_dict, output='residual'):
"""
"""
# first update instrument from input parameters
full_params = np.asarray(self.full_params)
full_params[self.flags] = reduced_params

self.instr.update_from_parameter_list(full_params[:self.npi])
self.params = full_params[self.npi:]

# need this for dsp
bmatx = self.plane_data.latVecOps['B']
wlen = self.instr.beam_wavelength

# working with Patrick's pick dicts
pick_angs_dict = data_dict['picks']
pick_xys_dict = data_dict['pick_xys']
# tvec_c = np.asarray(
# data_dict['options']['tvec'], dtype=float
# ).flatten()

# build residual
retval = np.array([], dtype=float)
for det_key, panel in self.instr.detectors.items():
# !!! now grabbing this from picks
hkls_ref = np.asarray(data_dict['hkls'][det_key], dtype=int)
if len(hkls_ref) == 0:
continue

gvecs = np.dot(hkls_ref, bmatx.T)
dsp_ref = 1./xfcapi.rowNorm(gvecs)

pick_angs = pick_angs_dict[det_key]
pick_xys = pick_xys_dict[det_key]
assert len(pick_angs) == len(dsp_ref), "picks are wrong length"
assert len(pick_angs) == len(pick_xys), "pick xy data inconsistent"
# the data structure is:
# [x, y, tth, eta, h, k, l, dsp0]
# FIXME: clean this up!
pdata = [] # [xy_meas, tth_eta_meas, hkl]
for ir, hkld in enumerate(zip(hkls_ref, dsp_ref)):
npts = len(pick_angs[ir])
if npts > 0:
tth_eta_meas = np.atleast_2d(np.radians(pick_angs[ir]))
xy_meas = np.atleast_2d(pick_xys[ir])
pdata.append(
np.hstack(
[xy_meas,
tth_eta_meas,
np.tile(np.hstack(hkld), (npts, 1))]
)
)
if len(pdata) == 0:
raise RuntimeWarning("detector '%s' has no feature points"
% det_key)
continue
else:
# recast as array
pdata = np.vstack(pdata)

"""
Here is the strategy:
1. remap the feature points from raw cartesian to
(tth, eta) under the current mapping
2. use the lattice and hkls to calculate the ideal tth0
3. push the (tth0, eta) values back through the mapping to
raw cartesian coordinates
4. build residual on the measured and recalculated (x, y)
"""
# the data structure is:
# [x, y, tth, eta, h, k, l, dsp0]
#
# push measured (x, y) ring points through current mapping
# to (tth, eta)
meas_xy = pdata[:, :2]
updates_angles, _ = panel.cart_to_angles(
meas_xy,
tvec_s=self.instr.tvec,
apply_distortion=True
)

# derive ideal tth positions from additional ring point info
hkls = pdata[:, 4:7]
gvecs = np.dot(hkls, bmatx.T)
dsp0 = 1./np.sqrt(np.sum(gvecs*gvecs, axis=1))
tth0 = 2.*np.arcsin(0.5*wlen/dsp0)

if output == 'residual':
retval = np.append(
retval,
updates_angles[:, 0].flatten() - tth0.flatten()
)
# retval = np.append(
# retval,
# meas_xy.flatten() - calc_xy.flatten()
# )
elif output == 'model':
# !!! get eta from mapped markers rather than ref
# eta0 = pdata[:, -1]
eta0 = updates_angles[:, 1]

# map updated (tth0, eta0) back to cartesian coordinates
tth_eta = np.vstack([tth0, eta0]).T
calc_xy = panel.angles_to_cart(
tth_eta,
tvec_s=self.instr.tvec,
apply_distortion=True
)
retval = np.append(
retval,
calc_xy.flatten()
)
else:
raise RuntimeError("unrecognized output flag '%s'"
% output)

return retval

def residual(self, reduced_params, data_dict):
return self._evaluate(reduced_params, data_dict)

def model(self, reduced_params, data_dict):
return self._evaluate(reduced_params, data_dict, output='model')


class CompositeCalibration(object):
def __init__(self, instr, processed_picks):
def __init__(self, instr, processed_picks, img_dict):
self.instr = instr
self.original_instr = copy.deepcopy(instr)
self.npi = len(self.instr.calibration_parameters)
self.data = processed_picks
calibrator_list = []
Expand All @@ -689,7 +499,7 @@ def __init__(self, instr, processed_picks):
)
param_flags.append(lpflags)
calib = PowderCalibrator(
self.instr, pick_data['plane_data'], flags
self.instr, pick_data['plane_data'], img_dict, flags
)
params.append(calib.full_params[-calib.npe:])
calibrator_list.append(calib)
Expand Down Expand Up @@ -727,10 +537,59 @@ def reduced_params(self):
return self.full_params[self.flags]

def residual(self, reduced_params, pick_data_list):
# first update full parameter list
self.full_params[self.flags] = reduced_params
instr_params = self.full_params[:self.npi]
addtl_params = self.full_params[self.npi:]
# first update a copy of the full parameter list
full_params = np.array(self.full_params)
full_params[self.flags] = reduced_params
instr_params = full_params[:self.npi]
addtl_params = full_params[self.npi:]

def powder_residual(calib, these_reduced_params, data_dict):
# Convert our data_dict into the input data format that
# the powder calibrator is expecting.
calibration_data = {}
for det_key in data_dict['hkls']:
# MUST use original instrument to convert these coordinates
# so that we are consistent. This is because the instrument
# gets modified with each call to `residual()`.
panel = self.original_instr.detectors[det_key]

picks_list = []
for i, hkl in enumerate(data_dict['hkls'][det_key]):
picks = data_dict['picks'][det_key][i]
if not picks:
continue

# Each row consists of 8 columns:
# First two are measured x, y
# Third is unknown (appears unused in the PowderCalibrator)
# Four - six are the hkl indices
# Seven and Eight are unknown (appears unused here)

# We are going to convert our data into rows of 6 columns.
# We will fill in zero for the third column, and ignore
# the last two columns.

# Convert the angles to Cartesian
cartesian_picks = angles_to_cart(picks, panel)
repeated_zero = np.repeat([[0]], len(cartesian_picks),
axis=0)
repeated_hkl = np.repeat([hkl], len(cartesian_picks),
axis=0)
formatted = np.hstack((cartesian_picks, repeated_zero,
repeated_hkl))
picks_list.append(formatted)

calibration_data[det_key] = picks_list

return calib.residual(these_reduced_params, calibration_data)

def laue_residual(calib, these_reduced_params, data_dict):
return calib.residual(these_reduced_params, data_dict)

residual_funcs = {
PowderCalibrator: powder_residual,
LaueCalibrator: laue_residual,
}

# loop calibrators and collect residuals
ii = 0
Expand All @@ -744,12 +603,10 @@ def residual(self, reduced_params, pick_data_list):
# pull out reduced list
these_reduced_params = these_full_params[calib.flags]

# call to calibrator residual api with porper index into pick data
# call to calibrator residual api with proper index into pick data
f = residual_funcs[type(calib)]
residual.append(
calib.residual(
these_reduced_params,
pick_data_list[ical]
)
f(calib, these_reduced_params, pick_data_list[ical])
)

# advance alibrator extra parameter offset
Expand All @@ -759,18 +616,26 @@ def residual(self, reduced_params, pick_data_list):
return np.hstack(residual)


def run_calibration(picks, instr, materials):
def run_calibration(picks, instr, img_dict, materials):
enrich_pick_data(picks, instr, materials)

# Run composite calibration
instr_calibrator = CompositeCalibration(instr, picks)
instr_calibrator = CompositeCalibration(instr, picks, img_dict)

x0_comp = instr_calibrator.reduced_params()

# Compute resd0, as hexrd does
resd0 = instr_calibrator.residual(x0_comp, picks) # noqa: F841

x1, cox_x, infodict, mesg, ierr = leastsq(
instr_calibrator.residual, x0_comp, args=(picks, ),
factor=0.1, full_output=True
)
full_output=True
)

# Evaluate new residual
# This will update the parameters
resd1 = instr_calibrator.residual(x1, picks) # noqa: F841

return instr_calibrator


Expand Down

0 comments on commit b6289b3

Please sign in to comment.