From 11a7b450c22a66610be3630e9a6e4f9bfc10b88d Mon Sep 17 00:00:00 2001 From: Patrick Avery Date: Wed, 13 Sep 2023 16:12:47 -0500 Subject: [PATCH] The calibration classes were transferred to HEXRD They better belong in HEXRD, not here. The files are being transferred to HEXRD in hexrd/hexrd#557. Depends on: hexrd/hexrd#557 Signed-off-by: Patrick Avery --- hexrd/ui/calibration/calibration_runner.py | 9 +- hexrd/ui/calibration/calibrationutil.py | 542 ----------------- .../ui/calibration/pick_based_calibration.py | 569 +----------------- 3 files changed, 5 insertions(+), 1115 deletions(-) delete mode 100644 hexrd/ui/calibration/calibrationutil.py diff --git a/hexrd/ui/calibration/calibration_runner.py b/hexrd/ui/calibration/calibration_runner.py index 98b8557a2..25bb4884d 100644 --- a/hexrd/ui/calibration/calibration_runner.py +++ b/hexrd/ui/calibration/calibration_runner.py @@ -9,7 +9,9 @@ from PySide2.QtWidgets import QFileDialog, QMessageBox from hexrd.crystallography import hklToStr -from hexrd.fitting.calibration import InstrumentCalibrator, PowderCalibrator +from hexrd.fitting.calibration import ( + InstrumentCalibrator, LaueCalibrator, PowderCalibrator +) from hexrd.instrument import unwrap_h5_to_dict from hexrd.ui.calibration.auto import ( @@ -17,10 +19,7 @@ save_picks_to_overlay, ) from hexrd.ui.calibration.laue_auto_picker_dialog import LaueAutoPickerDialog -from hexrd.ui.calibration.pick_based_calibration import ( - LaueCalibrator, - run_calibration, -) +from hexrd.ui.calibration.pick_based_calibration import run_calibration from hexrd.ui.calibration.hkl_picks_tree_view_dialog import ( generate_picks_results, overlays_to_tree_format, HKLPicksTreeViewDialog, picks_cartesian_to_angles, tree_format_to_picks, diff --git a/hexrd/ui/calibration/calibrationutil.py b/hexrd/ui/calibration/calibrationutil.py deleted file mode 100644 index fb99fa6c3..000000000 --- a/hexrd/ui/calibration/calibrationutil.py +++ /dev/null @@ -1,542 +0,0 @@ -#! /usr/bin/env python -# ============================================================================= -# Copyright (c) 2012, Lawrence Livermore National Security, LLC. -# Produced at the Lawrence Livermore National Laboratory. -# Written by Joel Bernier and others. -# LLNL-CODE-529294. -# All rights reserved. -# -# This file is part of HEXRD. For details on dowloading the source, -# see the file COPYING. -# -# Please also see the file LICENSE. -# -# This program is free software; you can redistribute it and/or modify it under -# the terms of the GNU Lesser General Public License (as published by the Free -# Software Foundation) version 2.1 dated February 1999. -# -# This program is distributed in the hope that it will be useful, but -# WITHOUT ANY WARRANTY; without even the IMPLIED WARRANTY OF MERCHANTABILITY -# or FITNESS FOR A PARTICULAR PURPOSE. See the terms and conditions of the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU Lesser General Public -# License along with this program (see file LICENSE); if not, write to -# the Free Software Foundation, Inc., 59 Temple Place, Suite 330, -# Boston, MA 02111-1307 USA or visit . -# ============================================================================= - -""" -Author: J. V. Bernier -Date: 16 Jan 2019 -""" - -# ============================================================================= -# # IMPORTS -# ============================================================================= - -import pickle as cpl - -import numpy as np - -import h5py - -import yaml - -# import ipywidgets as widgets -# from IPython.display import display - -from skimage import io -from skimage.exposure import rescale_intensity -from skimage.filters.edges import binary_erosion -from skimage.restoration import denoise_bilateral - -from hexrd import constants as cnst -from hexrd import instrument -from hexrd import matrixutil as mutil -from hexrd import rotations as rot -from hexrd import material - -from scipy.ndimage.morphology import binary_fill_holes -from scipy.optimize import least_squares, leastsq - - -# %% -# ============================================================================= -# # Parameters and local funcs -# ============================================================================= - -r2d = 180. / np.pi -d2r = np.pi / 180. -pi = np.pi -piby2 = 0.5*np.pi - -sigma_to_FWHM = 2*np.sqrt(2*np.log(2)) - -__reflInfo_dtype = [ - ('iRefl', int), - ('hkl', (int, 3)), - ('intensity', (float, 2)), - ('energy', float), - ('predAngles', (float, 2)), - ('measAngles', (float, 2)), - ('measXY', (float, 2)), - ] - - -class DummyCrystal(object): - def __init__(self, tth_list, delth_tth=np.radians(5.)): - self._tth = np.array(tth_list) - self._tThWidth = delth_tth - - @property - def tth(self): - return self._tth - - @tth.setter - def tth(self, x): - self._tth = np.radians(x) - - @property - def tThWidth(self): - return self._tThWidth - - @tThWidth.setter - def tThWidth(self, x): - self._tThWidth = x - - def getTTh(self): - return self.tth - - def getTThRanges(self): - tth_lo = self.getTTh() - 0.5*self.tThWidth - tth_hi = self.getTTh() + 0.5*self.tThWidth - return np.vstack([tth_lo, tth_hi]).T - - def getMergedRanges(self, cullDupl=False): - """ - return indices and ranges for specified planeData, merging where - there is overlap based on the tThWidth and line positions - """ - tThs = self.getTTh() - tThRanges = self.getTThRanges() - - # if you end exlcusions in a doublet (or multiple close rings) - # then this will 'fail'. May need to revisit... - nonoverlap_nexts = np.hstack( - (tThRanges[:-1, 1] < tThRanges[1:, 0], True) - ) - iHKLLists = [] - mergedRanges = [] - hklsCur = [] - tThLoIdx = 0 - tThHiCur = 0. - for iHKL, nonoverlapNext in enumerate(nonoverlap_nexts): - tThHi = tThRanges[iHKL, -1] - if not nonoverlapNext: - tth_diff = abs(tThs[iHKL] - tThs[iHKL + 1]) - if cullDupl and tth_diff < cnst.sqrt_epsf: - continue - else: - hklsCur.append(iHKL) - tThHiCur = tThHi - else: - hklsCur.append(iHKL) - tThHiCur = tThHi - iHKLLists.append(hklsCur) - mergedRanges.append([tThRanges[tThLoIdx, 0], tThHiCur]) - tThLoIdx = iHKL + 1 - hklsCur = [] - return iHKLLists, mergedRanges - - -# mask setting -def det_panel_mask(instr, img_dict, tolerance=1e-6): - """ - use small values surrounding image plate to set panel buffers - """ - for key, panel in instr.detectors.items(): - img = img_dict[key] - bimg = binary_fill_holes(img > tolerance) - mask = binary_erosion(bimg, iterations=3) - panel.panel_buffer = mask - - -def load_images(img_stem, ip_keys, - threshold=None, - denoise=False, - normalize=False): - img_dict = dict.fromkeys(ip_keys) - for ip_key in ip_keys: - this_img = io.imread(img_stem % ip_key.upper()) - if threshold is not None: - this_img[this_img < threshold] = 0. - if denoise: - this_img = np.array( - rescale_intensity( - denoise_bilateral(this_img, - multichannel=False, - sigma_spatial=1.1, - bins=2**16), - out_range=np.uint16), - dtype=np.uint16 - ) - if normalize: - this_img = rescale_intensity(this_img, out_range=(-1., 1.)) - img_dict[ip_key] = this_img - return img_dict - - -def log_scale_img(img): - img = np.array(img, dtype=float) - np.min(img) + 1. - return np.log(img) - - -# Material instantiation -# FIXME: these two functions are out of date! -def make_matl(mat_name, sgnum, lparms, hkl_ssq_max=500): - matl = material.Material(mat_name) - matl.sgnum = sgnum - matl.latticeParameters = lparms - matl.hklMax = hkl_ssq_max - - nhkls = len(matl.planeData.exclusions) - matl.planeData.set_exclusions(np.zeros(nhkls, dtype=bool)) - return matl - - -# crystallography data extraction from cPickle arhive -def load_plane_data(cpkl, key): - with open(cpkl, 'rb') as matf: - mat_list = cpl.load(matf) - pd = dict(zip([i.name for i in mat_list], mat_list))[key].planeData - pd.exclusions = np.zeros_like(pd.exclusions, dtype=bool) - return pd - - -# Tilt utilities -def convert_tilt(zxz_angles): - tilt = np.radians(zxz_angles) - rmat = rot.make_rmat_euler(tilt, 'zxz', extrinsic=False) - phi, n = rot.angleAxisOfRotMat(rmat) - return phi*n.flatten() - - -# pareser for simulation results -def parse_laue_simulation(sim_dict): - """ - !!!: assumes for single grain - ???: could eventually add another loop... - """ - gid = 0 - - # output dictionaries for each IP - valid_xy = dict.fromkeys(sim_dict) - valid_hkls = dict.fromkeys(sim_dict) - valid_energy = dict.fromkeys(sim_dict) - valid_angs = dict.fromkeys(sim_dict) - for ip_key, sim_results in sim_dict.items(): - # expand results for convenience - xy_det, hkls_in, angles, dspacing, energy = sim_results - valid_xy[ip_key] = [] - valid_hkls[ip_key] = [] - valid_energy[ip_key] = [] - valid_angs[ip_key] = [] - for gid in range(len(xy_det)): - # find valid reflections - valid_refl = ~np.isnan(xy_det[gid][:, 0]) - - valid_xy_tmp = xy_det[gid][valid_refl, :] - - # cull duplicates - dupl = mutil.findDuplicateVectors(valid_xy_tmp.T, tol=1e-4) - - # find hkls and angles to feed patchs - valid_xy[ip_key].append(valid_xy_tmp[dupl[1], :]) - valid_hkls[ip_key].append(hkls_in[gid][:, valid_refl][:, dupl[1]]) - valid_energy[ip_key].append(energy[gid][valid_refl]) - valid_angs[ip_key].append(angles[gid][valid_refl, :][dupl[1], :]) - - """ - # !!! not working for now - # need xy coords and pixel sizes - if distortion is not None: - valid_xy = distortion[0](valid_xy, - distortion[1], - invert=True) - """ - return valid_xy, valid_hkls, valid_energy, valid_angs - - -# Objective function for Laue fitting -def sxcal_obj_func(plist_fit, plist_full, param_flags, - instr, meas_xy, hkls_idx, - bmat, energy_cutoffs, - sim_only=False, - return_value_flag=None): - """ - Objective function for Laue-based fitting. - - - energy_cutoffs are [minEnergy, maxEnergy] where min/maxEnergy can be lists - - """ - npi_tot = len(instr.calibration_parameters) - - # fill out full parameter list - # !!! no scaling for now - plist_full[param_flags] = plist_fit - - plist_instr = plist_full[:npi_tot] - grain_params = [plist_full[npi_tot:], ] - - # update instrument - instr.update_from_parameter_list(plist_instr) - - # beam vector - bvec = instr.beam_vector - - # right now just stuck on the end and assumed - # to all be the same length... FIX THIS - calc_xy = {} - calc_ang = {} - npts_tot = 0 - for det_key, panel in instr.detectors.items(): - # counter - npts_tot += len(meas_xy[det_key]) - - # Simulate Laue pattern: - # returns xy_det, hkls_in, angles, dspacing, energy - sim_results = panel.simulate_laue_pattern( - [hkls_idx[det_key], bmat], - minEnergy=energy_cutoffs[0], maxEnergy=energy_cutoffs[1], - grain_params=grain_params, - beam_vec=bvec - ) - - calc_xy_tmp = sim_results[0][0] - calc_angs_tmp = sim_results[2][0] - - idx = ~np.isnan(calc_xy_tmp[:, 0]) - calc_xy[det_key] = calc_xy_tmp[idx, :] - calc_ang[det_key] = calc_angs_tmp[idx, :] - pass - - # return values - if sim_only: - retval = {} - for det_key in calc_xy.keys(): - # ??? calc_xy is always 2-d - retval[det_key] = [calc_xy[det_key], calc_ang[det_key]] - else: - meas_xy_all = [] - calc_xy_all = [] - for det_key in meas_xy.keys(): - meas_xy_all.append(meas_xy[det_key]) - calc_xy_all.append(calc_xy[det_key]) - pass - meas_xy_all = np.vstack(meas_xy_all) - calc_xy_all = np.vstack(calc_xy_all) - - diff_vecs_xy = calc_xy_all - meas_xy_all - retval = diff_vecs_xy.flatten() - if return_value_flag == 1: - retval = sum(abs(retval)) - elif return_value_flag == 2: - denom = npts_tot - len(plist_fit) - 1. - if denom != 0: - nu_fac = 1. / denom - else: - nu_fac = 1. - nu_fac = 1 / (npts_tot - len(plist_fit) - 1.) - retval = nu_fac * sum(retval**2) - return retval - - -# Calibration function -def calibrate_instrument_from_laue( - instr, grain_params, meas_xy, bmat, hkls_idx, - energy_cutoffs, param_flags=None, - xtol=cnst.sqrt_epsf, ftol=cnst.sqrt_epsf, - factor=1., sim_only=False, use_robust_lsq=False): - """ - """ - npi = len(instr.calibration_parameters) - - pnames = [ - '{:>24s}'.format('beam_energy'), - '{:>24s}'.format('beam_azimuth'), - '{:>24s}'.format('beam_polar_angle'), - '{:>24s}'.format('chi'), - '{:>24s}'.format('tvec_s[0]'), - '{:>24s}'.format('tvec_s[1]'), - '{:>24s}'.format('tvec_s[2]'), - ] - - for det_key, panel in instr.detectors.items(): - pnames += [ - '{:>24s}'.format('%s tilt[0]' % det_key), - '{:>24s}'.format('%s tilt[1]' % det_key), - '{:>24s}'.format('%s tilt[2]' % det_key), - '{:>24s}'.format('%s tvec[0]' % det_key), - '{:>24s}'.format('%s tvec[1]' % det_key), - '{:>24s}'.format('%s tvec[2]' % det_key), - ] - if panel.distortion is not None: - # FIXME: hard-coded distortion kludge - for dp in panel.distortion[1]: - pnames += ['{:>24s}'.format('%s distortion[0]' % det_key), ] - - pnames += [ - '{:>24s}'.format('crystal tilt[0]'), - '{:>24s}'.format('crystal tilt[1]'), - '{:>24s}'.format('crystal tilt[2]'), - '{:>24s}'.format('crystal tvec[0]'), - '{:>24s}'.format('crystal tvec[1]'), - '{:>24s}'.format('crystal tvec[2]'), - '{:>24s}'.format('crystal vinv[0]'), - '{:>24s}'.format('crystal vinv[1]'), - '{:>24s}'.format('crystal vinv[2]'), - '{:>24s}'.format('crystal vinv[3]'), - '{:>24s}'.format('crystal vinv[4]'), - '{:>24s}'.format('crystal vinv[5]'), - ] - - # reset parameter flags for instrument as specified - if param_flags is None: - param_flags_full = instr.calibration_flags - param_flags = np.hstack( - [param_flags_full, - np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=bool)] - ) - else: - # will throw an AssertionError if wrong length - assert(len(param_flags) == npi + 12) - instr.calibration_flags = param_flags[:npi] - - # set tilt mapping to ZXZ - # FIXME: input parameter? - # rme = rot.RotMatEuler(np.zeros(3), 'zxz', extrinsic=False) - # instr.tilt_calibration_mapping = rme - - # munge energy cutoffs - if hasattr(energy_cutoffs[0], '__len__'): - energy_cutoffs[0] = [0.5*i for i in energy_cutoffs[0]] - energy_cutoffs[1] = [1.5*i for i in energy_cutoffs[1]] - else: - energy_cutoffs[0] = 0.5*energy_cutoffs[0] - energy_cutoffs[1] = 1.5*energy_cutoffs[1] - - # grab relevant parameters - # will yield: - # 0 beam wavelength - # 1 beam azimuth - # 2 beam polar angle - # 3 chi - # 4:7 tvec_s - # panel_0 tilt, tvec, - # panel_1 tilt, tvec, - # ... - # panel_n tilt, tvec, - # grain_parameters - plist_i = instr.calibration_parameters - plist_full = np.hstack([plist_i, grain_params]) - plist_fit = plist_full[param_flags] - fit_args = (plist_full, param_flags, - instr, meas_xy, hkls_idx, - bmat, energy_cutoffs) - if sim_only: - return sxcal_obj_func( - plist_fit, plist_full, param_flags, - instr, meas_xy, hkls_idx, - bmat, energy_cutoffs, - sim_only=True) - else: - print("Set up to refine:") - for i in np.where(param_flags)[0]: - print("\t%s = %1.7e" % (pnames[i], plist_full[i])) - resd = sxcal_obj_func( - plist_fit, plist_full, param_flags, - instr, meas_xy, hkls_idx, - bmat, energy_cutoffs) - print("Initial SSR: %f" % (np.sqrt(np.sum(resd*resd)))) - - # run optimization - if use_robust_lsq: - result = least_squares( - sxcal_obj_func, plist_fit, args=fit_args, - xtol=xtol, ftol=ftol, - loss='soft_l1', method='trf' - ) - x = result.x - resd = result.fun - mesg = result.message - ierr = result.status - else: - # do least squares problem - x, cov_x, infodict, mesg, ierr = leastsq( - sxcal_obj_func, plist_fit, args=fit_args, - factor=factor, xtol=xtol, ftol=ftol, - full_output=1 - ) - resd = infodict['fvec'] - if ierr not in [1, 2, 3, 4]: - raise RuntimeError("solution not found: ierr = %d" % ierr) - else: - print("INFO: optimization fininshed successfully with ierr=%d" - % ierr) - print("INFO: %s" % mesg) - - # ??? output message handling? - fit_params = plist_full - fit_params[param_flags] = x - - print("Final parameter values:") - for i in np.where(param_flags)[0]: - print("\t%s = %1.7e" % (pnames[i], fit_params[i])) - print("Final SSR: %f" % (np.sqrt(np.sum(resd*resd)))) - # run simulation with optimized results - sim_final = sxcal_obj_func( - x, plist_full, param_flags, - instr, meas_xy, hkls_idx, - bmat, energy_cutoffs, - sim_only=True) - - ''' - # ??? reset instrument here? - instr.beam_vector = instrument.calc_beam_vec( - fit_params[0], fit_params[1]) - ii = npi # offset to where the panel parameters start - for det_key, panel in instr.detectors.items(): - panel.tilt = convert_tilt(fit_params[ii:ii + 3]) - panel.tvec = fit_params[ii + 3:ii + 6] - ii += npp - pass - ''' - return fit_params, resd, sim_final - - -# peak fitting -def gaussian_1d(p, x): - func = p[0]*np.exp(-(x-p[1])**2/2/p[2]**2) + p[3] - return func - - -def gaussian_2d(p, data): - shape = data.shape - x, y = np.meshgrid(range(shape[1]), range(shape[0])) - func = p[0]*np.exp( - -(p[1]*(x-p[4])*(x-p[4]) - + p[2]*(x-p[4])*(y-p[5]) - + p[3]*(y-p[5])*(y-p[5])) - ) + p[6]*(x-p[4]) + p[7]*(y-p[5]) + p[8] - return func.flatten() - data.flatten() - - -def gaussian_2d_int(y, x, *p): - func = p[0]*np.exp( - -(p[1]*(x-p[4])*(x-p[4]) - + p[2]*(x-p[4])*(y-p[5]) - + p[3]*(y-p[5])*(y-p[5])) - ) - return func.flatten() diff --git a/hexrd/ui/calibration/pick_based_calibration.py b/hexrd/ui/calibration/pick_based_calibration.py index 92271cf2d..562ab309e 100644 --- a/hexrd/ui/calibration/pick_based_calibration.py +++ b/hexrd/ui/calibration/pick_based_calibration.py @@ -1,30 +1,7 @@ -# -*- coding: utf-8 -*- -""" -Created on Tue Sep 29 14:20:48 2020 - -@author: berni -""" - -import copy - import numpy as np - -from scipy import ndimage -from scipy.integrate import nquad from scipy.optimize import leastsq -from skimage import filters -from skimage.feature import blob_log - -from hexrd.ui.calibration.calibrationutil import ( - 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 +from hexrd.fitting.calibration import CompositeCalibration def enrich_pick_data(picks, instr, materials): @@ -79,550 +56,6 @@ def enrich_pick_data(picks, instr, materials): pick_data[data_key][det_key] = pdl -# %% CLASSES - -class LaueCalibrator(object): - calibrator_type = 'laue' - _nparams = 12 - - def __init__(self, instr, plane_data, grain_params, flags, - min_energy=5., max_energy=25.): - self._instr = instr - self._plane_data = copy.deepcopy(plane_data) - self._plane_data.wavelength = self._instr.beam_energy # force - self._params = np.asarray(grain_params, dtype=float).flatten() - assert len(self._params) == self._nparams, \ - "grain parameters must have %d elements" % self._nparams - self._full_params = np.hstack( - [self._instr.calibration_parameters, self._params] - ) - assert len(flags) == len(self._full_params), \ - "flags must have %d elements; you gave %d" \ - % (len(self._full_params), len(flags)) - self._flags = flags - self._energy_cutoffs = [min_energy, max_energy] - - @property - def instr(self): - return self._instr - - @property - def plane_data(self): - self._plane_data.wavelength = self.energy_cutoffs[-1] - self._plane_data.exclusions = None - 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.params): - raise RuntimeError("params must have %d elements" - % len(self.params)) - self._params = 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] - - @property - def energy_cutoffs(self): - return self._energy_cutoffs - - @energy_cutoffs.setter - def energy_cutoffs(self, x): - assert len(x) == 2, "input must have 2 elements" - assert x[1] > x[0], "first element must be < than second" - self._energy_cutoffs = x - - def _autopick_points(self, raw_img_dict, tth_tol=5., eta_tol=5., - npdiv=2, do_smoothing=True, smoothing_sigma=2, - use_blob_detection=True, blob_threshold=0.25, - fit_peaks=True, min_peak_int=1., fit_tth_tol=0.1): - """ - - - Parameters - ---------- - raw_img_dict : TYPE - DESCRIPTION. - tth_tol : TYPE, optional - DESCRIPTION. The default is 5.. - eta_tol : TYPE, optional - DESCRIPTION. The default is 5.. - npdiv : TYPE, optional - DESCRIPTION. The default is 2. - do_smoothing : TYPE, optional - DESCRIPTION. The default is True. - smoothing_sigma : TYPE, optional - DESCRIPTION. The default is 2. - use_blob_detection : TYPE, optional - DESCRIPTION. The default is True. - blob_threshold : TYPE, optional - DESCRIPTION. The default is 0.25. - fit_peaks : TYPE, optional - DESCRIPTION. The default is True. - - Returns - ------- - None. - - """ - labelStructure = ndimage.generate_binary_structure(2, 1) - rmat_s = np.eye(3) # !!! forcing to identity - omega = 0. # !!! same ^^^ - - rmat_c = xfcapi.makeRotMatOfExpMap(self.params[:3]) - tvec_c = self.params[3:6] - # vinv_s = self.params[6:12] # !!!: patches don't take this yet - - # run simulation - # ???: could we get this from overlays? - laue_sim = self.instr.simulate_laue_pattern( - self.plane_data, - minEnergy=self.energy_cutoffs[0], - maxEnergy=self.energy_cutoffs[1], - rmat_s=None, grain_params=np.atleast_2d(self.params), - ) - - # loop over detectors for results - refl_dict = dict.fromkeys(self.instr.detectors) - for det_key, det in self.instr.detectors.items(): - det_config = det.config_dict( - chi=self.instr.chi, - tvec=self.instr.tvec, - beam_vector=self.instr.beam_vector - ) - - xy_det, hkls, angles, dspacing, energy = laue_sim[det_key] - ''' - valid_xy = [] - valid_hkls = [] - valid_angs = [] - valid_energy = [] - ''' - # !!! not necessary to loop over grains since we can only handle 1 - # for gid in range(len(xy_det)): - gid = 0 - # find valid reflections - valid_refl = ~np.isnan(xy_det[gid][:, 0]) - valid_xy = xy_det[gid][valid_refl, :] - valid_hkls = hkls[gid][:, valid_refl] - valid_angs = angles[gid][valid_refl, :] - valid_energy = energy[gid][valid_refl] - # pass - - # make patches - refl_patches = xrdutil.make_reflection_patches( - det_config, - valid_angs, det.angularPixelSize(valid_xy), - rmat_c=rmat_c, tvec_c=tvec_c, - tth_tol=tth_tol, eta_tol=eta_tol, - npdiv=npdiv, quiet=True) - - reflInfoList = [] - img = raw_img_dict[det_key] - native_area = det.pixel_area - num_patches = len(valid_angs) - meas_xy = np.nan*np.ones((num_patches, 2)) - meas_angs = np.nan*np.ones((num_patches, 2)) - for iRefl, patch in enumerate(refl_patches): - # check for overrun - irow = patch[-1][0] - jcol = patch[-1][1] - if np.any([irow < 0, irow >= det.rows, - jcol < 0, jcol >= det.cols]): - continue - if not np.all( - det.clip_to_panel( - np.vstack([patch[1][0].flatten(), - patch[1][1].flatten()]).T - )[1] - ): - continue - # use nearest interpolation - spot_data = img[irow, jcol] * patch[3] * npdiv**2 / native_area - spot_data -= np.amin(spot_data) - patch_size = spot_data.shape - - sigmax = 0.25*np.min(spot_data.shape) * fwhm_to_sigma - - # optional gaussian smoothing - if do_smoothing: - spot_data = filters.gaussian(spot_data, smoothing_sigma) - - if use_blob_detection: - spot_data_scl = 2.*spot_data/np.max(spot_data) - 1. - - # Compute radii in the 3rd column. - blobs_log = blob_log(spot_data_scl, - min_sigma=2, - max_sigma=min(sigmax, 20), - num_sigma=10, - threshold=blob_threshold, - overlap=0.1) - numPeaks = len(blobs_log) - else: - labels, numPeaks = ndimage.label( - spot_data > np.percentile(spot_data, 99), - structure=labelStructure - ) - slabels = np.arange(1, numPeaks + 1) - tth_edges = patch[0][0][0, :] - eta_edges = patch[0][1][:, 0] - delta_tth = tth_edges[1] - tth_edges[0] - delta_eta = eta_edges[1] - eta_edges[0] - if numPeaks > 0: - peakId = iRefl - if use_blob_detection: - coms = blobs_log[:, :2] - else: - coms = np.array( - ndimage.center_of_mass( - spot_data, labels=labels, index=slabels - ) - ) - if numPeaks > 1: - # - center = np.r_[spot_data.shape]*0.5 - com_diff = coms - np.tile(center, (numPeaks, 1)) - closest_peak_idx = np.argmin( - np.sum(com_diff**2, axis=1) - ) - # - else: - closest_peak_idx = 0 - pass # end multipeak conditional - # - coms = coms[closest_peak_idx] - # - if fit_peaks: - sigm = 0.2*np.min(spot_data.shape) - if use_blob_detection: - sigm = min(blobs_log[closest_peak_idx, 2], sigm) - y0, x0 = coms.flatten() - ampl = float(spot_data[int(y0), int(x0)]) - # y0, x0 = 0.5*np.array(spot_data.shape) - # ampl = np.max(spot_data) - a_par = c_par = 0.5/float(sigm**2) - b_par = 0. - bgx = bgy = 0. - bkg = np.min(spot_data) - params = [ampl, - a_par, b_par, c_par, - x0, y0, bgx, bgy, bkg] - # - result = leastsq(gaussian_2d, params, args=(spot_data)) - # - fit_par = result[0] - # - coms = np.array([fit_par[5], fit_par[4]]) - ''' - print("%s, %d, (%.2f, %.2f), (%d, %d)" - % (det_key, iRefl, coms[0], coms[1], - patch_size[0], patch_size[1])) - ''' - row_cen = fit_tth_tol * patch_size[0] - col_cen = fit_tth_tol * patch_size[1] - if np.any( - [coms[0] < row_cen, - coms[0] >= patch_size[0] - row_cen, - coms[1] < col_cen, - coms[1] >= patch_size[1] - col_cen] - ): - continue - if (fit_par[0] < min_peak_int): - continue - - # intensities - spot_intensity, int_err = nquad( - gaussian_2d_int, - [[0., 2.*y0], [0., 2.*x0]], - args=fit_par) - pass - com_angs = np.hstack([ - tth_edges[0] + (0.5 + coms[1])*delta_tth, - eta_edges[0] + (0.5 + coms[0])*delta_eta - ]) - - # grab intensities - if not fit_peaks: - if use_blob_detection: - spot_intensity = 10 - max_intensity = 10 - else: - spot_intensity = np.sum( - spot_data[labels == slabels[closest_peak_idx]] - ) - max_intensity = np.max( - spot_data[labels == slabels[closest_peak_idx]] - ) - else: - max_intensity = np.max(spot_data) - # need xy coords - # !!! forcing ome = 0. -- could be inconsistent with rmat_s - cmv = np.atleast_2d(np.hstack([com_angs, omega])) - gvec_c = xfcapi.anglesToGVec( - cmv, - chi=self.instr.chi, - rMat_c=rmat_c, - bHat_l=self.instr.beam_vector) - new_xy = xfcapi.gvecToDetectorXY( - gvec_c, - det.rmat, rmat_s, rmat_c, - det.tvec, self.instr.tvec, tvec_c, - beamVec=self.instr.beam_vector) - meas_xy[iRefl, :] = new_xy - if det.distortion is not None: - meas_xy[iRefl, :] = det.distortion.apply_inverse( - meas_xy[iRefl, :] - ) - meas_angs[iRefl, :] = com_angs - else: - peakId = -999 - # - spot_intensity = np.nan - max_intensity = np.nan - pass - reflInfoList.append([peakId, valid_hkls[:, iRefl], - (spot_intensity, max_intensity), - valid_energy[iRefl], - valid_angs[iRefl, :], - meas_angs[iRefl, :], - meas_xy[iRefl, :]]) - pass - reflInfo = np.array( - [tuple(i) for i in reflInfoList], - dtype=reflInfo_dtype) - refl_dict[det_key] = reflInfo - - # !!! ok, here is where we would populated the data_dict from refl_dict - return refl_dict - - def _evaluate(self, reduced_params, data_dict): - """ - """ - # 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:] - - # grab reflection data from picks input - pick_hkls_dict = dict.fromkeys(self.instr.detectors) - pick_xys_dict = dict.fromkeys(self.instr.detectors) - for det_key in self.instr.detectors: - # find valid reflections and recast hkls to int - xys = data_dict['pick_xys'][det_key] - hkls = np.asarray(data_dict['hkls'][det_key], dtype=int) - - valid_idx = ~np.isnan(xys[:, 0]) - - # fill local dicts - pick_hkls_dict[det_key] = np.atleast_2d(hkls[valid_idx, :]).T - pick_xys_dict[det_key] = np.atleast_2d(xys[valid_idx, :]) - - return pick_hkls_dict, pick_xys_dict - - def residual(self, reduced_params, data_dict): - # need this for laue obj - bmatx = self.plane_data.latVecOps['B'] - pick_hkls_dict, pick_xys_dict = self._evaluate( - reduced_params, data_dict - ) - # munge energy cutoffs - energy_cutoffs = np.r_[0.5, 1.5] * np.asarray(self.energy_cutoffs) - - return sxcal_obj_func( - reduced_params, self.full_params, self.flags, - self.instr, pick_xys_dict, pick_hkls_dict, - bmatx, energy_cutoffs - ) - - def model(self, reduced_params, data_dict): - # need this for laue obj - bmatx = self.plane_data.latVecOps['B'] - pick_hkls_dict, pick_xys_dict = self._evaluate( - reduced_params, data_dict, - ) - - return sxcal_obj_func( - reduced_params, self.full_params, self.flags, - self.instr, pick_xys_dict, pick_hkls_dict, - bmatx, self.energy_cutoffs, - sim_only=True - ) - - -class CompositeCalibration(object): - 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 = [] - params = [] - param_flags = [] - for pick_data in processed_picks: - if pick_data['type'] == 'powder': - # flags for calibrator - lpflags = [i[1] for i in pick_data['refinements']] - flags = np.hstack( - [self.instr.calibration_flags, lpflags] - ) - param_flags.append(lpflags) - - kwargs = { - 'instr': self.instr, - 'plane_data': pick_data['plane_data'], - 'img_dict': img_dict, - 'flags': flags, - 'tth_distortion': pick_data['tth_distortion'], - } - calib = PowderCalibrator(**kwargs) - - params.append(calib.full_params[-calib.npe:]) - calibrator_list.append(calib) - - elif pick_data['type'] == 'laue': - # flags for calibrator - gparams = pick_data['options']['crystal_params'] - min_energy = pick_data['options']['min_energy'] - max_energy = pick_data['options']['max_energy'] - - gpflags = [i[1] for i in pick_data['refinements']] - flags = np.hstack( - [self.instr.calibration_flags, gpflags] - ) - param_flags.append(gpflags) - calib = LaueCalibrator( - self.instr, pick_data['plane_data'], - gparams, flags, - min_energy=min_energy, max_energy=max_energy - ) - params.append(calib.full_params[-calib.npe:]) - calibrator_list.append(calib) - - self.calibrators = calibrator_list - self.params = np.hstack(params) - self.param_flags = np.hstack(param_flags) - self.full_params = np.hstack( - [self.instr.calibration_parameters, self.params] - ) - self.flags = np.hstack( - [self.instr.calibration_flags, self.param_flags] - ) - - def reduced_params(self): - return self.full_params[self.flags] - - def residual(self, reduced_params, pick_data_list): - # 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 - residual = [] - for ical, calib in enumerate(self.calibrators): - # make copy offull params for this calibrator - these_full_params = np.hstack( - [instr_params, addtl_params[ii:ii + calib.npe]] - ) - - # pull out reduced list - these_reduced_params = these_full_params[calib.flags] - - # call to calibrator residual api with proper index into pick data - f = residual_funcs[type(calib)] - residual.append( - f(calib, these_reduced_params, pick_data_list[ical]) - ) - - # advance alibrator extra parameter offset - ii += calib.npe - - # return single hstacked residual - return np.hstack(residual) - - def run_calibration(picks, instr, img_dict, materials): enrich_pick_data(picks, instr, materials)