From c3039fc26626ad748127316a38a7afd26fd8c618 Mon Sep 17 00:00:00 2001 From: Kevin Lewis Date: Mon, 29 Jul 2024 12:06:54 -0400 Subject: [PATCH] Remove scripts folder --- .../detector_templates}/Hydra_Feb19.yml | 0 scripts/calibrate_from_powder.py | 591 ------------- scripts/calibrate_from_powder.py.bak | 691 --------------- scripts/calibrate_from_rotation_series.py | 612 ------------- scripts/calibration_mockup.py | 469 ---------- scripts/convert_instrument_config.py | 161 ---- scripts/pinhole_2theta.py | 836 ------------------ scripts/plot_instrument_schematic.py | 179 ---- scripts/powder_calibration.py | 249 ------ scripts/preprocess_tomo_mask.py | 214 ----- scripts/process_nf.py | 364 -------- scripts/tiffs_from_h5.py | 64 -- 12 files changed, 4430 deletions(-) rename {scripts => hexrd/resources/detector_templates}/Hydra_Feb19.yml (100%) delete mode 100644 scripts/calibrate_from_powder.py delete mode 100644 scripts/calibrate_from_powder.py.bak delete mode 100644 scripts/calibrate_from_rotation_series.py delete mode 100644 scripts/calibration_mockup.py delete mode 100644 scripts/convert_instrument_config.py delete mode 100644 scripts/pinhole_2theta.py delete mode 100644 scripts/plot_instrument_schematic.py delete mode 100644 scripts/powder_calibration.py delete mode 100644 scripts/preprocess_tomo_mask.py delete mode 100644 scripts/process_nf.py delete mode 100644 scripts/tiffs_from_h5.py diff --git a/scripts/Hydra_Feb19.yml b/hexrd/resources/detector_templates/Hydra_Feb19.yml similarity index 100% rename from scripts/Hydra_Feb19.yml rename to hexrd/resources/detector_templates/Hydra_Feb19.yml diff --git a/scripts/calibrate_from_powder.py b/scripts/calibrate_from_powder.py deleted file mode 100644 index 6a94a47e0..000000000 --- a/scripts/calibrate_from_powder.py +++ /dev/null @@ -1,591 +0,0 @@ -#!/usr/bin/env python2 -# -*- coding: utf-8 -*- -""" -Created on Thu Sep 14 12:04:31 2017 - -@author: bernier2 -""" -import os - -import glob - -import matplotlib.gridspec as gridspec -from matplotlib import pyplot as plt - -import numpy as np - -from scipy.optimize import leastsq, least_squares -try: - import dill as cpl -except(ImportError): - import pickle as cpl - -from skimage.exposure import equalize_adapthist, rescale_intensity -from skimage.morphology import disk, binary_erosion - -import yaml - -from hexrd import imageseries -from hexrd import instrument -from hexrd import material -from hexrd.fitting import fitpeak -from hexrd.rotations import \ - angleAxisOfRotMat, angles_from_rmat_xyz, make_rmat_euler -from hexrd.transforms.xfcapi import mapAngle -from hexrd.xrdutil import make_reflection_patches - - -# plane data -def load_pdata(cpkl, key): - with open(cpkl, "rb") as matf: - mat_list = cpl.load(matf) - return dict(list(zip([i.name for i in mat_list], mat_list)))[key].planeData - - -# images -def load_images(yml): - return imageseries.open(yml, format="image-files") - - -# instrument -def load_instrument(yml): - with open(yml, 'r') as f: - icfg = yaml.safe_load(f) - return instrument.HEDMInstrument(instrument_config=icfg) - - -# build a material on the fly -def make_matl(mat_name, sgnum, lparms, hkl_ssq_max=50): - 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 - - -# %% -# ============================================================================= -# USER INPUT -# ============================================================================= - -# dirs -working_dir = '/Users/Shared/APS/PUP_AFRL_Feb19' -image_dir = os.path.join(working_dir, 'image_data') - -samp_name = 'ceria_cal' -scan_number = 0 - -raw_data_dir_template = os.path.join( - image_dir, - 'raw_images_%s_%06d-%s.yml' -) - -instrument_filename = 'Hydra_Feb19.yml' - -#mat_filename = 'materials.hexrd' -#mat_key = 'ruby' -#plane_data = load_pdata(os.path.join(working_dir, mat_filename), mat_key) -#plane_data.lparms = np.r_[2.508753] -matl = make_matl('ceria', 225, [5.411102, ]) - -# tolerances for patches -tth_tol = 0.15 -eta_tol = 5. - -tth_max = 9. - -# peak fit type -#pktype = 'gaussian' -pktype = 'pvoigt' - -# Try it -use_robust_optimization = False - -#refinement_type = "translation_only" -refinement_type = "translation_and_tilt" -#refinement_type = "beam_and_translation" - - -# %% -# ============================================================================= -# IMAGESERIES -# ============================================================================= - -# for applying processing options (flips, dark subtretc...) -PIS = imageseries.process.ProcessedImageSeries - -# load instrument -#instr = load_instrument(os.path.join(working_dir, instrument_filename)) -instr = load_instrument(instrument_filename) -det_keys = list(instr.detectors.keys()) - -# hijack panel buffer -edge_buff = 5 -for k, v in list(instr.detectors.items()): - try: - mask = np.load(mask_template % k) - except(NameError): - mask = np.ones((v.rows, v.cols), dtype=bool) - mask[:edge_buff, :] = False - mask[:, :edge_buff] = False - mask[-edge_buff:, :] = False - mask[:, -edge_buff:] = False - v.panel_buffer = mask - - -# grab imageseries filenames -file_names = glob.glob( - raw_data_dir_template % (samp_name, scan_number, '*') -) -check_files_exist = [os.path.exists(file_name) for file_name in file_names] -if not np.all(check_files_exist): - raise RuntimeError("files don't exist!") - -img_dict = dict.fromkeys(det_keys) -for k, v in img_dict.items(): - file_name = glob.glob( - raw_data_dir_template % (samp_name, scan_number, k) - ) - ims = load_images(file_name[0]) - img_dict[k] = ims[0] - -# plane data munging -matl.beamEnergy = instr.beam_energy -plane_data = matl.planeData -if tth_tol is not None: - plane_data.tThWidth = np.radians(tth_tol) - -if tth_max is not None: - plane_data.exclusions = None - plane_data.tThMax = np.radians(tth_max) - -# %% -# ============================================================================= -# INSTRUMENT -# ============================================================================= - - -# get powder line profiles -# -# output is as follows: -# -# patch_data = {:[ringset_index][patch_index]} -# patch_data[0] = [two_theta_edges, ref_eta] -# patch_data[1] = [intensities] -# -powder_lines = instr.extract_line_positions( - plane_data, img_dict, - tth_tol=np.degrees(plane_data.tThWidth), eta_tol=eta_tol, - npdiv=2, - collapse_eta=False, collapse_tth=False, - do_interpolation=True) - - -# ideal tth -tth_ideal = plane_data.getTTh() -tth0 = [] -for idx in plane_data.getMergedRanges()[0]: - tth0.append(tth_ideal[idx[0]]) - -# GRAND LOOP OVER PATCHES -rhs = dict.fromkeys(det_keys) -for det_key in det_keys: - rhs[det_key] = [] - panel = instr.detectors[det_key] - for i_ring, ringset in enumerate(powder_lines[det_key]): - tmp = [] - for angs, intensities in ringset: - tth_centers = np.average( - np.vstack([angs[0][:-1], angs[0][1:]]), - axis=0) - eta_ref = angs[1] - int1d = np.sum(np.array(intensities).squeeze(), axis=0) - """ - DARREN: FIT [tth_centers, intensities[0]] HERE - - RETURN TTH0 - rhs.append([tth0, eta_ref]) - """ - p0 = fitpeak.estimate_pk_parms_1d( - tth_centers, int1d, pktype - ) - - p = fitpeak.fit_pk_parms_1d( - p0, tth_centers, int1d, pktype - ) - # !!! this is where we can kick out bunk fits - tth_ref = plane_data.getTTh()[i_ring] - tth_meas = p[1] - if p[0] < 0.1 or abs(tth_meas - tth_ref) > np.radians(tth_tol/2.): - continue - xy_meas = panel.angles_to_cart([[tth_meas, eta_ref], ]) - tmp.append( - np.hstack([xy_meas.squeeze(), tth_meas, tth0[i_ring], eta_ref]) - ) - rhs[det_key].append(np.vstack(tmp)) - rhs[det_key] = np.array(rhs[det_key]) - -# %% plot fit poistions -fig, ax = plt.subplots(2, 2) -fig_row, fig_col = np.unravel_index(np.arange(instr.num_panels), (2, 2)) - -ifig = 0 -for det_key, panel in instr.detectors.items(): - all_pts = np.vstack(rhs[det_key]) - ''' - pimg = equalize_adapthist( - rescale_intensity(img_dict[det_key], out_range=(0, 1)), - 10, clip_limit=0.2) - ''' - pimg = np.array(img_dict[det_key], dtype=float) - pimg[~panel.panel_buffer] = np.nan - ax[fig_row[ifig], fig_col[ifig]].imshow( - pimg, - vmin=np.percentile(img_dict[det_key], 5), - vmax=np.percentile(img_dict[det_key], 90), - cmap=plt.cm.bone_r - ) - ideal_angs, ideal_xys = panel.make_powder_rings(plane_data, - delta_eta=eta_tol) - rijs = panel.cartToPixel(np.vstack(ideal_xys)) - ax[fig_row[ifig], fig_col[ifig]].plot(rijs[:, 1], rijs[:, 0], 'cx') - ax[fig_row[ifig], fig_col[ifig]].set_title(det_key) - rijs = panel.cartToPixel(all_pts[:, :2]) - ax[fig_row[ifig], fig_col[ifig]].plot(rijs[:, 1], rijs[:, 0], 'm+') - ax[fig_row[ifig], fig_col[ifig]].set_title(det_key) - ifig += 1 - -# %% -if refinement_type == "translation_only": - x0 = [] - for k, v in list(instr.detectors.items()): - x0.append(v.tvec) - x0 = np.hstack(x0) - - def multipanel_powder_objfunc(param_list, data_dict, instr): - """ - """ - npp = 3 - ii = 0 - jj = ii + npp - resd = [] - for det_key, panel in list(instr.detectors.items()): - # strip params for this panel - params = param_list[ii:jj] - - # translation - panel.tvec = np.asarray(params).flatten() - - # advance counters - ii += npp - jj += npp - - pdata = np.vstack(data_dict[det_key]) - if len(pdata) > 0: - - calc_xy = panel.angles_to_cart(pdata[:, -2:]) - - resd.append( - (pdata[:, :2].flatten() - calc_xy.flatten()) - ) - else: - continue - return np.hstack(resd) -elif refinement_type == "translation_and_tilt": - # build parameter list - x0 = [] - for k, v in instr.detectors.items(): - tilt = np.degrees(angles_from_rmat_xyz(v.rmat)) - print("XYZ tilt angles for '%s': (%.2f, %.2f, %.2f)" - % tuple([k, *tilt])) - x0.append(np.hstack([tilt[:2], v.tvec])) - x0 = np.hstack(x0) - - def multipanel_powder_objfunc(param_list, data_dict, instr): - """ - """ - npp = 5 - ii = 0 - jj = ii + npp - resd = [] - ''' - resd = 0 - ''' - for det_key, panel in instr.detectors.items(): - # strip params for this panel - params = param_list[ii:jj] - - # tilt first - tilt = np.degrees(angles_from_rmat_xyz(panel.rmat)) - tilt[:2] = params[:2] - rmat = make_rmat_euler( - np.radians(tilt), 'xyz', extrinsic=True - ) - phi, n = angleAxisOfRotMat(rmat) - panel.tilt = phi*n.flatten() - - # translation - panel.tvec = params[2:] - - # advance counters - ii += npp - jj += npp - - ''' - for rdata in data_dict[det_key]: - cangs, cgvecs = panel.cart_to_angles(rdata[:, :2]) - resd += np.sum(np.abs(cangs[:, 0] - np.mean(cangs[:, 0]))) - ''' - pdata = np.vstack(data_dict[det_key]) - if len(pdata) > 0: - calc_xy = panel.angles_to_cart(pdata[:, -2:]) - - resd.append( - (pdata[:, :2].flatten() - calc_xy.flatten()) - ) - else: - continue - return np.hstack(resd) -elif refinement_type == "beam_and_translation": - # build parameter list - x0 = list(instrument.calc_angles_from_beam_vec(instr.beam_vector)) - for k, v in list(instr.detectors.items()): - x0.append(v.tvec) - x0 = np.hstack(x0) - - def multipanel_powder_objfunc(param_list, data_dict, instr): - """ - """ - # update beam vector - instr.beam_vector = param_list[:2] - - # number of parameters per panel - npp = 3 - - # initialize counters and residual - ii = 2 - jj = ii + npp - resd = [] - #resd = 0 - for det_key, panel in list(instr.detectors.items()): - # strip params for this panel - params = param_list[ii:jj] - - # update translation - panel.tvec = params - - # advance counters - ii += npp - jj += npp - - #for rdata in data_dict[det_key]: - # cangs, cgvecs = panel.cart_to_angles(rdata[:, :2]) - # resd += np.sum(np.abs(cangs[:, 0] - np.mean(cangs[:, 0]))) - pdata = np.vstack(data_dict[det_key]) - if len(pdata) > 0: - calc_xy = panel.angles_to_cart(pdata[:, -2:]) - - - resd.append( - (pdata[:, :2].flatten() - calc_xy.flatten()) - ) - else: - continue - return np.hstack(resd) # resd - - -# %% - -det_key = 'ge2' - -fig, ax = plt.subplots() - -rid = 0 -ii = 0 -for ring_data in powder_lines[det_key][rid]: - tth_edges = ring_data[0][0] - tth_centers = 0.5*np.sum(np.vstack([tth_edges[:-1], tth_edges[1:]]), axis=0) - eta = np.degrees(mapAngle(ring_data[0][1], (0, 2*np.pi))) - intensities = np.sum(np.array(ring_data[1]).squeeze(), axis=0) - if ii < 1: - iprev = 0 - else: - iprev = 0.25*np.max(np.sum(powder_lines[det_key][rid][ii - 1][1], axis=0)) - ax.plot(tth_centers, intensities + iprev) - ii += 1 - -p0 = fitpeak.estimate_pk_parms_1d( - tth_centers, intensities, pktype - ) - -p = fitpeak.fit_pk_parms_1d( - p0, tth_centers, intensities, pktype - ) - -if pktype == 'pvoigt': - fp0 = fitpeak.pkfuncs.pvoigt1d(p0, tth_centers) - fp1 = fitpeak.pkfuncs.pvoigt1d(p, tth_centers) -elif pktype == 'gaussian': - fp0 = fitpeak.pkfuncs.gaussian1d(p0, tth_centers) - fp1 = fitpeak.pkfuncs.gaussian1d(p, tth_centers) - -# %% - -initial_guess = np.array(x0) - -if use_robust_optimization: - oresult = least_squares( - multipanel_powder_objfunc, x0, args=(rhs, instr), - method='trf', loss='soft_l1' - ) - x1 = oresult['x'] -else: - x1, cox_x, infodict, mesg, ierr = leastsq( - multipanel_powder_objfunc, x0, args=(rhs, instr), - full_output=True - ) -resd0 = multipanel_powder_objfunc(x0, rhs, instr) -resd1 = multipanel_powder_objfunc(x1, rhs, instr) - -delta_r = sum(resd0**2)/float(len(resd0)) - sum(resd1**2)/float(len(resd1)) - -if delta_r > 0: - print(('OPTIMIZATION SUCCESSFUL\nfinal ssr: %f' % sum(resd1**2))) - print(('delta_r: %f' % delta_r)) - instr.write_config(instrument_filename) - x0 = np.array(x1) -else: - print('no improvement in residual!!!') - -# %% - - -# ============================================================================= -# %% PLOTTING -# ============================================================================= - -cmap = plt.cm.bone_r - -#img = np.log(img_dict[det_key] + (1. - np.min(img_dict[det_key]))) -img = equalize_adapthist(img_dict[det_key], 20, clip_limit=0.2) -panel = instr.detectors[det_key] - -pcfg = panel.config_dict(instr.chi, instr.tvec) -extent = [-0.5*panel.col_dim, 0.5*panel.col_dim, - -0.5*panel.row_dim, 0.5*panel.row_dim] -pr = panel.make_powder_rings(plane_data, - delta_tth=np.degrees(plane_data.tThWidth), - delta_eta=eta_tol) - -tth_eta = pr[0][rid] -xy_det = pr[1][rid] -rpatches = make_reflection_patches( - pcfg, tth_eta, panel.angularPixelSize(xy_det), - omega=None, tth_tol=np.degrees(plane_data.tThWidth), eta_tol=eta_tol, - distortion=panel.distortion, - npdiv=2, - beamVec=instr.beam_vector) - -# Create 2x2 sub plots -widths = [1, 3] -heights = [1, 3] -fig = plt.figure(constrained_layout=True) -gs = gridspec.GridSpec(ncols=2, nrows=2, - figure=fig, - width_ratios=widths, - height_ratios=heights) -ax0 = fig.add_subplot(gs[0, 0]) # row 0, col 0 -ax1 = fig.add_subplot(gs[0, 1]) # row 0, col 1 -ax2 = fig.add_subplot(gs[1, :]) # row 1, span all columns - -#vmin = np.percentile(img, 50) -#vmax = np.percentile(img, 99) -vmin = None -vmax = None -img[~panel.panel_buffer] = np.nan - -ax2.imshow(img, - vmin=vmin, - vmax=vmax, - cmap=cmap, - extent=extent) - -for rp in rpatches: - ''' - pts = np.dot(panel.rmat, - np.vstack([rp[1][0].flatten(), - rp[1][1].flatten(), - np.zeros(rp[1][0].size)]) - ) - px = pts[0, :] - py = pts[1, :] - ''' - px = rp[1][0].flatten() - py = rp[1][1].flatten() - _, on_panel = panel.clip_to_panel(np.vstack([px, py]).T) - if np.any(~on_panel): - continue - else: - good_patch = [px, py] - ax2.plot(px, py, 'm.', markersize=0.1) -ax2.plot(good_patch[0], good_patch[1], 'c.', markersize=0.1) -aext = np.degrees( - [np.min(rp[0][0]), - np.max(rp[0][0]), - np.min(rp[0][1]), - np.max(rp[0][1])] -) -ax0.imshow(np.array(ring_data[1]).squeeze(), cmap=cmap, - extent=aext, aspect=.05, - vmin=np.percentile(ring_data[1], 50), - vmax=np.percentile(ring_data[1], 99)) -ax1.plot(np.degrees(tth_centers), intensities, 'k.-') -ax1.plot(np.degrees(tth_centers), fp1, 'r:') - -ax2.set_xlabel(r'detector $X$ [mm]') -ax2.set_ylabel(r'detector $Y$ [mm]') - -ax0.set_xlabel(r'patch $2\theta$') -ax0.set_ylabel(r'patch $\eta$') -ax0.axis('auto') - -ax1.set_xlabel(r'patch $2\theta$') -ax1.set_ylabel(r'summed intensity [arb]') - - -## ============================================================================= -## %% ERROR PLOTTING -## ============================================================================= -# -#tths = plane_data.getTTh() -#shkls = plane_data.getHKLs(asStr=True) -#ref_etas = np.radians(np.linspace(0, 360, num=181, endpoint=True)) -# -#ii = 1 -#for tth, shkl in zip(tths, shkls): -# #fig, ax = plt.subplots(1, 1, figsize=(12, 6)) -# #fig.suptitle('{%s}' % shkl) -# fig = plt.figure(ii) -# ax = fig.gca() -# ax.set_xlim(-180, 180) -# ax.set_ylim(-0.001, 0.001) -# ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%1.1e')) -# ax.set_xlabel(r'azimuth, $\eta$ [deg]', size=18) -# ax.set_ylabel(r'Relative error, $\frac{2\theta-2\theta_0}{2\theta_0}$', size=18) -# ax.grid(True) -# for det_key, rdata in rhs.iteritems(): -# if len(rdata) > 0: -# idx = rdata[:, 3] == tth -# if np.any(idx): -# dtth = (rdata[idx, 2] - tth*np.ones(sum(idx)))/tth -# ax.scatter(np.degrees(rdata[idx, -1]), dtth, c='b', s=12) -# ii += 1 -# -## %% -#for i, shkl in enumerate(shkls): -# fig = plt.figure(i+1) -# ax = fig.gca() -# plt.savefig("CeO2_tth_errors_%s.png" % '_'.join(shkl.split()), dpi=180) diff --git a/scripts/calibrate_from_powder.py.bak b/scripts/calibrate_from_powder.py.bak deleted file mode 100644 index 1bc9cc24b..000000000 --- a/scripts/calibrate_from_powder.py.bak +++ /dev/null @@ -1,691 +0,0 @@ -#!/usr/bin/env python2 -# -*- coding: utf-8 -*- -""" -Created on Thu Sep 14 12:04:31 2017 - -@author: bernier2 -""" -import os - -import glob - -import matplotlib.gridspec as gridspec -from matplotlib import pyplot as plt - -import numpy as np - -from scipy.optimize import leastsq, least_squares, minimize -try: - import dill as cpl -except(ImportError): - import pickle as cpl - -from skimage.exposure import equalize_adapthist, rescale_intensity -from skimage.morphology import disk, binary_erosion - -import yaml - -from hexrd import imageseries -from hexrd import instrument -from hexrd.fitting import fitpeak -from hexrd.rotations import \ - angleAxisOfRotMat, angles_from_rmat_zxz, make_rmat_euler -from hexrd.transforms.xfcapi import mapAngle, makeRotMatOfExpMap -from hexrd.xrdutil import make_reflection_patches - - -# plane data -def load_pdata(cpkl, key): - with file(cpkl, "r") as matf: - mat_list = cpl.load(matf) - return dict(list(zip([i.name for i in mat_list], mat_list)))[key].planeData - - -# images -def load_images(yml): - return imageseries.open(yml, format="image-files") - - -# instrument -def load_instrument(yml): - with file(yml, 'r') as f: - icfg = yaml.load(f) - return instrument.HEDMInstrument(instrument_config=icfg) - - -# %% -# ============================================================================= -# USER INPUT -# ============================================================================= - -#samp_name = 'N170209-003' -#samp_name = 'N170419-002' -samp_name = 'N190516-001' - -single_panel = False - -if single_panel: - raw_data_dir_template = './%s*.h5' - instrument_filename = 'pilatus6M_20190324.yml' -else: - raw_data_dir_template = './%s*images.yml' - #instrument_filename = 'tardis_N170209_BCC.yml' - #instrument_filename = 'tardis_N170419.yml' - instrument_filename = 'tardis_N190516-001.yml' - #instrument_filename = 'new-instrument-0.yml' - -load_processed = True # load hdf5 with all processed images in it - -mat_filename = './materials.hexrd' -#mat_key = 'BCC Mg 307GPa' -#mat_key = 'FCC Mg 550GPa' -#mat_key = 'SC Mg' -mat_key = 'FCC Al' - -mask_template = './%s_%%s_lmask.npy' % samp_name - -if load_processed: - pass -else: - # !!!: hard coded options for each GE for Oct 2017 - # CAVEAT: these keys MUST match instrument file! - # - # in case of dark file: - # - # panel_opts = [('dark', , ('flip', 'h'), ] - # - # note that panel options are applied in L-R order from the list. - panel_opts = None - -# tolerances for patches -tth_tol = 4. -eta_tol = 2. - -tth_max = 65. - -# peak fit type -pktype = 'gaussian' -# pktype = 'pvoigt' - -# Try it -use_robust_optimization = False - -#refinement_type = "translation_only" -refinement_type = "translation_and_tilt" -#refinement_type = "beam_and_translation" - - -# %% -# ============================================================================= -# IMAGESERIES -# ============================================================================= - -# for applying processing options (flips, dark subtretc...) -PIS = imageseries.process.ProcessedImageSeries - -# load instrument -instr = load_instrument(instrument_filename) -det_keys = list(instr.detectors.keys()) - -# hijack panel buffer -edge_buff = 10 -for k, v in instr.detectors.items(): - mask = np.load(mask_template % k) - mask[:edge_buff, :] = False - mask[:, :edge_buff] = False - mask[-edge_buff:, :] = False - mask[:, -edge_buff:] = False - v.panel_buffer = mask - - -# grab imageseries filenames -file_names = glob.glob( - raw_data_dir_template % samp_name -) -check_files_exist = [os.path.exists(file_name) for file_name in file_names] -if not np.all(check_files_exist): - raise RuntimeError("files don't exist!") - -img_dict = dict.fromkeys(det_keys) -if load_processed: - ims = load_images(file_names[0]) - img_dict = dict.fromkeys(ims.metadata['panels']) - for i_det, det_key in enumerate(ims.metadata['panels']): - img_dict[det_key] = ims[i_det] -else: - for file_name in file_names: - ims = load_images(file_name) - this_key = file_name.split('/')[-1].split('-')[1].upper() - if panel_opts[this_key] is not None: - img_dict[this_key] = PIS(ims, panel_opts[this_key])[0] - else: - img_dict[this_key] = ims[0] - - -# %% -# ============================================================================= -# INSTRUMENT -# ============================================================================= - -plane_data = load_pdata(mat_filename, mat_key) -#plane_data.lparms = np.r_[2.508753] -if tth_tol is not None: - plane_data.tThWidth = np.radians(tth_tol) - -if tth_max is not None: - plane_data.exclusions = None - plane_data.tThMax = np.radians(tth_max) - -# get powder line profiles -# -# output is as follows: -# -# patch_data = {:[ringset_index][patch_index]} -# patch_data[0] = [two_theta_edges, ref_eta] -# patch_data[1] = [intensities] -# -powder_lines = instr.extract_line_positions( - plane_data, img_dict, - tth_tol=np.degrees(plane_data.tThWidth), eta_tol=eta_tol, - npdiv=2, - collapse_eta=False, collapse_tth=False, - do_interpolation=True) - - -# ideal tth -tth_ideal = plane_data.getTTh() -tth0 = [] -for idx in plane_data.getMergedRanges()[0]: - tth0.append(tth_ideal[idx[0]]) - -# GRAND LOOP OVER PATCHES -rhs = dict.fromkeys(det_keys) -for det_key in det_keys: - rhs[det_key] = [] - panel = instr.detectors[det_key] - for i_ring, ringset in enumerate(powder_lines[det_key]): - tmp = [] - for angs, intensities in ringset: - tth_centers = np.average( - np.vstack([angs[0][:-1], angs[0][1:]]), - axis=0) - eta_ref = angs[1] - int1d = np.sum(np.array(intensities).squeeze(), axis=0) - """ - DARREN: FIT [tth_centers, intensities[0]] HERE - - RETURN TTH0 - rhs.append([tth0, eta_ref]) - """ - p0 = fitpeak.estimate_pk_parms_1d( - tth_centers, int1d, pktype - ) - - p = fitpeak.fit_pk_parms_1d( - p0, tth_centers, int1d, pktype - ) - # !!! this is where we can kick out bunk fits - tth_ref = plane_data.getTTh()[i_ring] - if p[0] < 0.1 or abs(p[1] - tth_ref) > np.radians(tth_tol/3.): - continue - tth_meas = p[1] - # - xy_meas = panel.angles_to_cart([[tth_meas, eta_ref], ]) - tmp.append( - np.hstack([xy_meas.squeeze(), tth_meas, tth0[i_ring], eta_ref]) - ) - rhs[det_key].append(np.vstack(tmp)) - rhs[det_key] = np.array(rhs[det_key]) - -# %% plot fit poistions -fig, ax = plt.subplots(2, 1) -ifig = 0 -for det_key, panel in instr.detectors.items(): - all_pts = np.vstack(rhs[det_key]) - pimg = equalize_adapthist( - rescale_intensity(img_dict[det_key], out_range=(0, 1)), - 10, clip_limit=0.2) - pimg[~panel.panel_buffer] = np.nan - ax[ifig].imshow( - pimg - ) - ideal_angs, ideal_xys = panel.make_powder_rings(plane_data, - delta_eta=eta_tol) - rijs = panel.cartToPixel(np.vstack(ideal_xys)) - ax[ifig].plot(rijs[:, 1], rijs[:, 0], 'c.') - ax[ifig].set_title(det_key) - rijs = panel.cartToPixel(all_pts[:, :2]) - ax[ifig].plot(rijs[:, 1], rijs[:, 0], 'mo') - ax[ifig].set_title(det_key) - ifig += 1 - -# %% -if refinement_type == "translation_only": - x0 = [] - for k, v in instr.detectors.items(): - x0.append(v.tvec) - x0 = np.hstack(x0) - - def multipanel_powder_objfunc(param_list, data_dict, instr): - """ - """ - npp = 3 - ii = 0 - jj = ii + npp - resd = [] - for det_key, panel in instr.detectors.items(): - # strip params for this panel - params = param_list[ii:jj] - - # translation - panel.tvec = np.asarray(params).flatten() - - # advance counters - ii += npp - jj += npp - - pdata = np.vstack(data_dict[det_key]) - if len(pdata) > 0: - - calc_xy = panel.angles_to_cart(pdata[:, -2:]) - - resd.append( - (pdata[:, :2].flatten() - calc_xy.flatten()) - ) - else: - continue - return np.hstack(resd) -elif refinement_type == "translation_and_tilt": - # build parameter list - x0 = [] - for k, v in instr.detectors.items(): - tilt = np.degrees(angles_from_rmat_zxz(v.rmat)) - x0.append(np.hstack([tilt[2], v.tvec])) - x0 = np.hstack(x0) - - def multipanel_powder_objfunc(param_list, data_dict, instr): - """ - """ - npp = 4 - ii = 0 - jj = ii + npp - resd = [] - ''' - resd = 0 - ''' - for det_key, panel in instr.detectors.items(): - # strip params for this panel - params = param_list[ii:jj] - - # tilt first - tilt = np.degrees(angles_from_rmat_zxz(panel.rmat)) - tilt[2] = params[0] - rmat = make_rmat_euler( - np.radians(tilt), 'zxz', extrinsic=False - ) - phi, n = angleAxisOfRotMat(rmat) - panel.tilt = phi*n - - # translation - panel.tvec = params[1:] - - # advance counters - ii += npp - jj += npp - - ''' - for rdata in data_dict[det_key]: - cangs, cgvecs = panel.cart_to_angles(rdata[:, :2]) - resd += np.sum(np.abs(cangs[:, 0] - np.mean(cangs[:, 0]))) - ''' - pdata = np.vstack(data_dict[det_key]) - if len(pdata) > 0: - calc_xy = panel.angles_to_cart(pdata[:, -2:]) - - resd.append( - (pdata[:, :2].flatten() - calc_xy.flatten()) - ) - else: - continue - return np.hstack(resd) -elif refinement_type == "beam_and_translation": - # build parameter list - x0 = list(instrument.calc_angles_from_beam_vec(instr.beam_vector)) - for k, v in instr.detectors.items(): - x0.append(v.tvec) - x0 = np.hstack(x0) - - def multipanel_powder_objfunc(param_list, data_dict, instr): - """ - """ - # update beam vector - instr.beam_vector = param_list[:2] - - # number of parameters per panel - npp = 3 - - # initialize counters and residual - ii = 2 - jj = ii + npp - resd = [] - #resd = 0 - for det_key, panel in instr.detectors.items(): - # strip params for this panel - params = param_list[ii:jj] - - # update translation - panel.tvec = params - - # advance counters - ii += npp - jj += npp - - #for rdata in data_dict[det_key]: - # cangs, cgvecs = panel.cart_to_angles(rdata[:, :2]) - # resd += np.sum(np.abs(cangs[:, 0] - np.mean(cangs[:, 0]))) - pdata = np.vstack(data_dict[det_key]) - if len(pdata) > 0: - calc_xy = panel.angles_to_cart(pdata[:, -2:]) - - - resd.append( - (pdata[:, :2].flatten() - calc_xy.flatten()) - ) - else: - continue - return np.hstack(resd) # resd - - -# %% - -if single_panel: - det_key = 'IP4' -else: - det_key = 'IP4' - -from hexrd.xrd.transforms import mapAngle - -fig, ax = plt.subplots() - -rid = 0 -ii = 0 -for ring_data in powder_lines[det_key][rid]: - tth_edges = ring_data[0][0] - tth_centers = 0.5*np.sum(np.vstack([tth_edges[:-1], tth_edges[1:]]), axis=0) - eta = np.degrees(mapAngle(ring_data[0][1], (0, 2*np.pi))) - intensities = np.sum(np.array(ring_data[1]).squeeze(), axis=0) - if ii < 1: - iprev = 0 - else: - iprev = 0.25*np.max(np.sum(powder_lines[det_key][rid][ii - 1][1], axis=0)) - ax.plot(tth_centers, intensities + iprev) - ii += 1 - -p0 = fitpeak.estimate_pk_parms_1d( - tth_centers, intensities, pktype - ) - -p = fitpeak.fit_pk_parms_1d( - p0, tth_centers, intensities, pktype - ) - -if pktype == 'pvoigt': - fp0 = fitpeak.pkfuncs.pvoigt1d(p0, tth_centers) - fp1 = fitpeak.pkfuncs.pvoigt1d(p, tth_centers) -elif pktype == 'gaussian': - fp0 = fitpeak.pkfuncs.gaussian1d(p0, tth_centers) - fp1 = fitpeak.pkfuncs.gaussian1d(p, tth_centers) - -# %% - -initial_guess = np.array(x0) - -if use_robust_optimization: - oresult = least_squares( - multipanel_powder_objfunc, x0, args=(rhs, instr), - method='trf', loss='soft_l1' - ) - x1 = oresult['x'] -else: - x1, cox_x, infodict, mesg, ierr = leastsq( - multipanel_powder_objfunc, x0, args=(rhs, instr), - full_output=True - ) -resd0 = multipanel_powder_objfunc(x0, rhs, instr) -resd1 = multipanel_powder_objfunc(x1, rhs, instr) - -delta_r = sum(resd0**2)/float(len(resd0)) - sum(resd1**2)/float(len(resd1)) - -if delta_r > 0: - print(('OPTIMIZATION SUCCESSFUL\nfinal ssr: %f' % sum(resd1**2))) - print(('delta_r: %f' % delta_r)) - instr.write_config(instrument_filename) - x0 = np.array(x1) -else: - print('no improvement in residual!!!') - -# %% - -delta_r = 1. -while delta_r > 0.1: - # get powder line profiles - # - # output is as follows: - # - # patch_data = {:[ringset_index][patch_index]} - # patch_data[0] = [two_theta_edges, ref_eta] - # patch_data[1] = [intensities] - # - powder_lines = instr.extract_line_positions( - plane_data, img_dict, - tth_tol=np.degrees(plane_data.tThWidth), eta_tol=eta_tol, - npdiv=2, - collapse_eta=False, collapse_tth=False, - do_interpolation=True) - - - # ideal tth - tth_ideal = plane_data.getTTh() - tth0 = [] - for idx in plane_data.getMergedRanges()[0]: - tth0.append(tth_ideal[idx[0]]) - - # GRAND LOOP OVER PATCHES - rhs = dict.fromkeys(det_keys) - for det_key in det_keys: - rhs[det_key] = [] - panel = instr.detectors[det_key] - for i_ring, ringset in enumerate(powder_lines[det_key]): - tmp = [] - for angs, intensities in ringset: - tth_centers = np.average( - np.vstack([angs[0][:-1], angs[0][1:]]), - axis=0) - eta_ref = angs[1] - int1d = np.sum(np.array(intensities).squeeze(), axis=0) - """ - DARREN: FIT [tth_centers, intensities[0]] HERE - - RETURN TTH0 - rhs.append([tth0, eta_ref]) - """ - p0 = fitpeak.estimate_pk_parms_1d( - tth_centers, int1d, pktype - ) - - p = fitpeak.fit_pk_parms_1d( - p0, tth_centers, int1d, pktype - ) - # !!! this is where we can kick out bunk fits - tth_ref = plane_data.getTTh()[i_ring] - if p[0] < 0.0001 or abs(p[1] - tth_ref) > np.radians(tth_tol/10.): - continue - tth_meas = p[1] - # - xy_meas = panel.angles_to_cart([[tth_meas, eta_ref], ]) - tmp.append( - np.hstack([xy_meas.squeeze(), tth_meas, tth0[i_ring], eta_ref]) - ) - rhs[det_key].append(np.vstack(tmp)) - - initial_guess = np.array(x0) - - if use_robust_optimization: - oresult = least_squares( - multipanel_powder_objfunc, x0, args=(rhs, instr), - method='trf', loss='soft_l1' - ) - x1 = oresult['x'] - else: - x1, cox_x, infodict, mesg, ierr = leastsq( - multipanel_powder_objfunc, x0, args=(rhs, instr), - full_output=True - ) - ''' - oresult = minimize(multipanel_powder_objfunc, x0, args=(rhs, instr), - method='BFGS', - ) - x1 = oresult['x'] - ''' - - resd0 = multipanel_powder_objfunc(x0, rhs, instr) - resd1 = multipanel_powder_objfunc(x1, rhs, instr) - - delta_r = resd0 - resd1 - #delta_r = sum(resd0**2)/float(len(resd0)) - sum(resd1**2)/float(len(resd1)) - - if delta_r > 0: - print(('OPTIMIZATION SUCCESSFUL\nfinal ssr: %f' % sum(resd1**2))) - print(('delta_r: %f' % delta_r)) - instr.write_config(instrument_filename) - x0 = np.array(x1) - else: - print('no improvement in residual!!!') - -# ============================================================================= -# %% PLOTTING -# ============================================================================= - -cmap = plt.cm.bone_r - -#img = np.log(img_dict[det_key] + (1. - np.min(img_dict[det_key]))) -img = equalize_adapthist(img_dict[det_key], 20, clip_limit=0.2) -panel = instr.detectors[det_key] - -pcfg = panel.config_dict(instr.chi, instr.tvec) -extent = [-0.5*panel.col_dim, 0.5*panel.col_dim, - -0.5*panel.row_dim, 0.5*panel.row_dim] -pr = panel.make_powder_rings(plane_data, - delta_tth=np.degrees(plane_data.tThWidth), - delta_eta=eta_tol) - -tth_eta = pr[0][rid] -xy_det = pr[1][rid] -rpatches = make_reflection_patches( - pcfg, tth_eta, panel.angularPixelSize(xy_det), - omega=None, tth_tol=np.degrees(plane_data.tThWidth), eta_tol=eta_tol, - distortion=panel.distortion, - npdiv=2, - beamVec=instr.beam_vector) - -# Create 2x2 sub plots -widths = [1, 3] -heights = [1, 3] -fig = plt.figure(constrained_layout=True) -gs = gridspec.GridSpec(ncols=2, nrows=2, - figure=fig, - width_ratios=widths, - height_ratios=heights) -ax0 = fig.add_subplot(gs[0, 0]) # row 0, col 0 -ax1 = fig.add_subplot(gs[0, 1]) # row 0, col 1 -ax2 = fig.add_subplot(gs[1, :]) # row 1, span all columns - -#vmin = np.percentile(img, 50) -#vmax = np.percentile(img, 99) -vmin = None -vmax = None -img[~panel.panel_buffer] = np.nan - -ax2.imshow(img, - vmin=vmin, - vmax=vmax, - cmap=cmap, - extent=extent) - -for rp in rpatches: - ''' - pts = np.dot(panel.rmat, - np.vstack([rp[1][0].flatten(), - rp[1][1].flatten(), - np.zeros(rp[1][0].size)]) - ) - px = pts[0, :] - py = pts[1, :] - ''' - px = rp[1][0].flatten() - py = rp[1][1].flatten() - _, on_panel = panel.clip_to_panel(np.vstack([px, py]).T) - if np.any(~on_panel): - continue - else: - good_patch = [px, py] - ax2.plot(px, py, 'm.', markersize=0.1) - -ax2.plot(good_patch[0], good_patch[1], 'c.', markersize=0.1) -aext = np.degrees( - [np.min(rp[0][0]), - np.max(rp[0][0]), - np.min(rp[0][1]), - np.max(rp[0][1])] -) -ax0.imshow(np.array(ring_data[1]).squeeze(), cmap=cmap, - extent=aext, aspect=.05, - vmin=np.percentile(ring_data[1], 50), - vmax=np.percentile(ring_data[1], 99)) -ax1.plot(np.degrees(tth_centers), intensities, 'k.-') -ax1.plot(np.degrees(tth_centers), fp1, 'r:') - -ax2.set_xlabel(r'detector $X$ [mm]') -ax2.set_ylabel(r'detector $Y$ [mm]') - -ax0.set_xlabel(r'patch $2\theta$') -ax0.set_ylabel(r'patch $\eta$') -ax0.axis('auto') - -ax1.set_xlabel(r'patch $2\theta$') -ax1.set_ylabel(r'summed intensity [arb]') - - -## ============================================================================= -## %% ERROR PLOTTING -## ============================================================================= -# -#tths = plane_data.getTTh() -#shkls = plane_data.getHKLs(asStr=True) -#ref_etas = np.radians(np.linspace(0, 360, num=181, endpoint=True)) -# -#ii = 1 -#for tth, shkl in zip(tths, shkls): -# #fig, ax = plt.subplots(1, 1, figsize=(12, 6)) -# #fig.suptitle('{%s}' % shkl) -# fig = plt.figure(ii) -# ax = fig.gca() -# ax.set_xlim(-180, 180) -# ax.set_ylim(-0.001, 0.001) -# ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%1.1e')) -# ax.set_xlabel(r'azimuth, $\eta$ [deg]', size=18) -# ax.set_ylabel(r'Relative error, $\frac{2\theta-2\theta_0}{2\theta_0}$', size=18) -# ax.grid(True) -# for det_key, rdata in rhs.iteritems(): -# if len(rdata) > 0: -# idx = rdata[:, 3] == tth -# if np.any(idx): -# dtth = (rdata[idx, 2] - tth*np.ones(sum(idx)))/tth -# ax.scatter(np.degrees(rdata[idx, -1]), dtth, c='b', s=12) -# ii += 1 -# -## %% -#for i, shkl in enumerate(shkls): -# fig = plt.figure(i+1) -# ax = fig.gca() -# plt.savefig("CeO2_tth_errors_%s.png" % '_'.join(shkl.split()), dpi=180) diff --git a/scripts/calibrate_from_rotation_series.py b/scripts/calibrate_from_rotation_series.py deleted file mode 100644 index 174233d9a..000000000 --- a/scripts/calibrate_from_rotation_series.py +++ /dev/null @@ -1,612 +0,0 @@ -ee -from __future__ import print_function - -import os - -import numpy as np - -try: - import dill as cpl -except(ImportError): - import cPickle as cpl - -import yaml - -from hexrd import config -from hexrd import constants as cnst -from hexrd import matrixutil as mutil -from hexrd import instrument -from hexrd.xrd import fitting -from hexrd.xrd import transforms_CAPI as xfcapi - -from scipy.optimize import leastsq, least_squares - -from matplotlib import pyplot as plt - - -# ============================================================================= -# %% *USER INPUT* -# ============================================================================= - -# hexrd yaml config file -cfg_filename = 'ruby_cal_config.yml' -block_id = 0 # only change this if you know what you are doing! - -# select which orientaion to use (in case of more than one...) -grain_id = 2 -override_grain_params = False - -# load previously saved exclusions -use_saved_exclusions = False -excl_filename = 'exclusion_index_%05d.npy' % grain_id - -overwrite_cfg = True - -# refit tolerances -n_pixels_tol = 2.0 -ome_tol = 2.0 - - -# ============================================================================= -# %% Local function definitions -# ============================================================================= - -# plane data -def load_plane_data(cpkl, key): - with file(cpkl, "r") as matf: - mat_list = cpl.load(matf) - return dict(zip([i.name for i in mat_list], mat_list))[key].planeData - - -# instrument -def load_instrument(yml): - with file(yml, 'r') as f: - icfg = yaml.load(f) - instr = instrument.HEDMInstrument(instrument_config=icfg) - for det_key in instr.detectors: - if 'saturation_level' in icfg['detectors'][det_key].keys(): - sat_level = icfg['detectors'][det_key]['saturation_level'] - print("INFO: Setting panel '%s' saturation level to %e" - % (det_key, sat_level)) - instr.detectors[det_key].saturation_level = sat_level - return instr - - -def calibrate_instrument_from_sx( - instr, grain_params, bmat, xyo_det, hkls_idx, - param_flags=None, dparam_flags=None, ome_period=None, - xtol=cnst.sqrt_epsf, ftol=cnst.sqrt_epsf, - factor=10., sim_only=False, use_robust_lsq=False): - """ - arguments xyo_det, hkls_idx are DICTs over panels - - !!! - distortion is still hosed... - Currently a dict of detector keys with - distortion[key] = [d_func, d_params, d_flags] - """ - pnames = [ - '{:>24s}'.format('wavelength'), - '{:>24s}'.format('chi'), - '{:>24s}'.format('tvec_s[0]'), - '{:>24s}'.format('tvec_s[1]'), - '{:>24s}'.format('tvec_s[2]'), - '{:>24s}'.format('expmap_c[0]'), - '{:>24s}'.format('expmap_c[0]'), - '{:>24s}'.format('expmap_c[0]'), - '{:>24s}'.format('tvec_c[0]'), - '{:>24s}'.format('tvec_c[1]'), - '{:>24s}'.format('tvec_c[2]'), - ] - - for det_key, panel in instr.detectors.iteritems(): - 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), - ] - - # now add distortion if - for det_key, panel in instr.detectors.iteritems(): - if panel.distortion is not None: - for j in range(len(panel.distortion[1])): - pnames.append( - '{:>24s}'.format('%s dparam[%d]' % (det_key, j)) - ) - - # reset parameter flags for instrument as specified - if param_flags is None: - param_flags = instr.param_flags - else: - # will throw an AssertionError if wrong length - instr.param_flags = param_flags - - det_keys = instr.detectors.keys() - - # re-map omegas if need be - if ome_period is not None: - for det_key in det_keys: - xyo_det[det_key][:, 2] = xfcapi.mapAngle( - xyo_det[det_key][:, 2], - ome_period - ) - # grain parameters - expmap_c = grain_params[:3] - tvec_c = grain_params[3:6] - vinv_s = grain_params[6:] - - plist_full = instr.calibration_params(expmap_c, tvec_c) - - dfuncs = {} - ndparams = {} - dparam_flags = [] - for det_key, panel in instr.detectors.iteritems(): - if panel.distortion is not None: - dfuncs[det_key] = panel.distortion[0] - # ititialize to all False... - ndparams[det_key] = len(panel.distortion[1]) - else: - dfuncs[det_key] = None - ndparams[det_key] = 0 - dparam_flags.append(np.zeros(ndparams[det_key], dtype=bool)) - dparam_flags = np.hstack(dparam_flags) - refine_flags = np.hstack([param_flags, dparam_flags]) - plist_fit = plist_full[refine_flags] - fit_args = (plist_full, param_flags, - dfuncs, dparam_flags, ndparams, - instr, xyo_det, hkls_idx, - bmat, vinv_s, ome_period, - instr.beam_vector, instr.eta_vector) - if sim_only: - return sxcal_obj_func( - plist_fit, plist_full, param_flags, - dfuncs, dparam_flags, ndparams, - instr, xyo_det, hkls_idx, - bmat, vinv_s, ome_period, - instr.beam_vector, instr.eta_vector, - sim_only=True) - else: - print("Set up to refine:") - for i in np.where(refine_flags)[0]: - print("\t%s = %1.7e" % (pnames[i], plist_full[i])) - - # 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[refine_flags] = x - - # run simulation with optimized results - sim_final = sxcal_obj_func( - x, plist_full, param_flags, - dfuncs, dparam_flags, ndparams, - instr, xyo_det, hkls_idx, - bmat, vinv_s, ome_period, - instr.beam_vector, instr.eta_vector, - sim_only=True) - - # ??? reset instrument here? - instr.beam_energy = cnst.keVToAngstrom(fit_params[0]) - instr.chi = fit_params[1] - instr.tvec = fit_params[2:5] - ii = 11 - for det_key, panel in instr.detectors.iteritems(): - panel.tilt = fit_params[ii:ii + 3] - panel.tvec = fit_params[ii + 3:ii + 6] - ii += 6 - # !!! use jj to do distortion? - - return fit_params, resd, sim_final - - -def sxcal_obj_func(plist_fit, plist_full, param_flags, - dfuncs, dparam_flags, ndparams, - instr, xyo_det, hkls_idx, - bmat, vinv_s, ome_period, - bvec, evec, - sim_only=False, return_value_flag=None): - """ - """ - # stack flags and force bool repr - refine_flags = np.array( - np.hstack([param_flags, dparam_flags]), - dtype=bool) - - # fill out full parameter list - # !!! no scaling for now - plist_full[refine_flags] = plist_fit - - # instrument quantities - wavelength = plist_full[0] - chi = plist_full[1] - tvec_s = plist_full[2:5] - - # calibration crystal quantities - rmat_c = xfcapi.makeRotMatOfExpMap(plist_full[5:8]) - tvec_c = plist_full[8:11] - - # right now just stuck on the end and assumed - # to all be the same length... FIX THIS - dparams_all = plist_full[-len(dparam_flags):] - xy_unwarped = {} - meas_omes = {} - calc_omes = {} - calc_xy = {} - ii = 11 # offset to start of panels... - jj = 0 - npts_tot = 0 - for det_key, panel in instr.detectors.iteritems(): - xy_unwarped[det_key] = xyo_det[det_key][:, :2] - npts_tot += len(xyo_det[det_key]) - dfunc = dfuncs[det_key] - len_these_dps = ndparams[det_key] - if dfunc is not None: # do unwarping - dparams = dparams_all[jj:jj + len_these_dps] - jj += len_these_dps - xy_unwarped[det_key] = dfunc(xy_unwarped[det_key], dparams) - meas_omes[det_key] = xyo_det[det_key][:, 2] - - # get these panel params for convenience - gparams = plist_full[ii:ii + 6] - - rmat_d = xfcapi.makeDetectorRotMat(gparams[:3]) - tvec_d = gparams[3:].reshape(3, 1) - - # transform G-vectors: - # 1) convert inv. stretch tensor from MV notation in to 3x3 - # 2) take reciprocal lattice vectors from CRYSTAL to SAMPLE frame - # 3) apply stretch tensor - # 4) normalize reciprocal lattice vectors in SAMPLE frame - # 5) transform unit reciprocal lattice vetors back to CRYSAL frame - gvec_c = np.dot(bmat, hkls_idx[det_key].T) - vmat_s = mutil.vecMVToSymm(vinv_s) - ghat_s = mutil.unitVector(np.dot(vmat_s, np.dot(rmat_c, gvec_c))) - ghat_c = np.dot(rmat_c.T, ghat_s) - - match_omes, calc_omes_tmp = fitting.matchOmegas( - xyo_det[det_key], hkls_idx[det_key].T, - chi, rmat_c, bmat, wavelength, - vInv=vinv_s, - beamVec=bvec, etaVec=evec, - omePeriod=ome_period) - - rmat_s_arr = xfcapi.make_oscill_rot_mat_array( - chi, np.ascontiguousarray(calc_omes_tmp) - ) - calc_xy_tmp = xfcapi.gvec_to_xy( - ghat_c.T, rmat_d, rmat_s_arr, rmat_c, - tvec_d, tvec_s, tvec_c - ) - if np.any(np.isnan(calc_xy_tmp)): - print("infeasible parameters: " - + "may want to scale back finite difference step size") - - calc_omes[det_key] = calc_omes_tmp - calc_xy[det_key] = calc_xy_tmp - - ii += 6 - - # return values - if sim_only: - retval = {} - for det_key in calc_xy.keys(): - # ??? calc_xy is always 2-d - retval[det_key] = np.vstack( - [calc_xy[det_key].T, calc_omes[det_key]] - ).T - else: - meas_xy_all = [] - calc_xy_all = [] - meas_omes_all = [] - calc_omes_all = [] - for det_key in xy_unwarped.keys(): - meas_xy_all.append(xy_unwarped[det_key]) - calc_xy_all.append(calc_xy[det_key]) - meas_omes_all.append(meas_omes[det_key]) - calc_omes_all.append(calc_omes[det_key]) - meas_xy_all = np.vstack(meas_xy_all) - calc_xy_all = np.vstack(calc_xy_all) - meas_omes_all = np.hstack(meas_omes_all) - calc_omes_all = np.hstack(calc_omes_all) - - diff_vecs_xy = calc_xy_all - meas_xy_all - diff_ome = xfcapi.angularDifference(calc_omes_all, meas_omes_all) - retval = np.hstack( - [diff_vecs_xy, - diff_ome.reshape(npts_tot, 1)] - ).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 - - -def parse_reflection_tables(cfg, instr, grain_id, refit_idx=None): - """ - make spot dictionaries - """ - hkls = {} - xyo_det = {} - idx_0 = {} - for det_key, panel in instr.detectors.iteritems(): - spots_filename = os.path.join( - cfg.analysis_dir, os.path.join( - det_key, 'spots_%05d.out' % grain_id - ) - ) - - # load pull_spots output table - gtable = np.loadtxt(spots_filename) - - # apply conditions for accepting valid data - valid_reflections = gtable[:, 0] >= 0 # is indexed - not_saturated = gtable[:, 6] < panel.saturation_level - print("INFO: panel '%s', grain %d" % (det_key, grain_id)) - print("INFO: %d of %d reflections are indexed" - % (sum(valid_reflections), len(gtable)) - ) - print("INFO: %d of %d" % (sum(not_saturated), sum(valid_reflections)) + - " valid reflections be are below" + - " saturation threshold of %d" % (panel.saturation_level) - ) - - # valid reflections index - if refit_idx is None: - idx = np.logical_and(valid_reflections, not_saturated) - idx_0[det_key] = idx - else: - idx = refit_idx[det_key] - idx_0[det_key] = idx - print("INFO: input reflection specify " + - "%d of %d total valid reflections" % (sum(idx), len(gtable)) - ) - - hkls[det_key] = gtable[idx, 2:5] - meas_omes = gtable[idx, 12].reshape(sum(idx), 1) - xyo_det[det_key] = np.hstack([gtable[idx, -2:], meas_omes]) - return hkls, xyo_det, idx_0 - - -# %% Initialization... - -# read config -cfg = config.open(cfg_filename)[block_id] - -# output for eta-ome maps as pickles -working_dir = cfg.working_dir -analysis_name = cfg.analysis_name -analysis_dir = cfg.analysis_dir -analysis_id = "%s_%s" % (analysis_name, cfg.material.active) - -# instrument -parfilename = cfg.instrument.parameters -instr = load_instrument(parfilename) -num_panels = instr.num_panels -det_keys = instr.detectors.keys() - -# plane data -plane_data = load_plane_data(cfg.material.definitions, cfg.material.active) -bmat = plane_data.latVecOps['B'] - -# the maximum pixel dimension in the instrument for plotting -max_pix_size = 0. -for panel in instr.detectors.itervalues(): - max_pix_size = max(max_pix_size, - max(panel.pixel_size_col, panel.pixel_size_col) - ) - -# grab omega period -# !!! data should be consistent -# !!! this is in degrees -ome_period = cfg.find_orientations.omega.period - -# load reflection tables from grain fit -hkls, xyo_det, idx_0 = parse_reflection_tables(cfg, instr, grain_id) - -# load grain parameters -if override_grain_params: - grain_parameters = np.loadtxt( - os.path.join(cfg.analysis_dir, 'grains.out'), - ndmin=2)[grain_id, 3:15] -else: - instr_cfg = yaml.load(open(cfg.instrument.parameters, 'r')) - try: - cal_crystal_dict = instr_cfg['calibration_crystal'] - if cal_crystal_dict['grain_id'] != grain_id: - raise RuntimeError("grain_id from instrument config is not %d" - % grain_id) - else: - grain_parameters = np.hstack( - [cal_crystal_dict['orientation'], - cal_crystal_dict['position'], - cal_crystal_dict['inv_stretch']] - ) - except(KeyError): - raise RuntimeError("instrument config lacks a calibration crystal") - - -# ============================================================================= -# %% plot initial guess -# ============================================================================= - -xyo_i = calibrate_instrument_from_sx( - instr, grain_parameters, bmat, xyo_det, hkls, - ome_period=np.radians(ome_period), sim_only=True -) - -for det_key, panel in instr.detectors.iteritems(): - fig, ax = plt.subplots(1, 2, sharex=True, sharey=False, figsize=(9, 5)) - fig.suptitle("detector %s" % det_key) - ax[0].plot(xyo_det[det_key][:, 0], xyo_det[det_key][:, 1], 'k.') - ax[0].plot(xyo_i[det_key][:, 0], xyo_i[det_key][:, 1], 'rx') - ax[0].grid(True) - ax[0].axis('equal') - ax[0].set_xlim(-0.5*panel.col_dim, 0.5*panel.col_dim) - ax[0].set_ylim(-0.5*panel.row_dim, 0.5*panel.row_dim) - ax[0].set_xlabel('detector X [mm]') - ax[0].set_ylabel('detector Y [mm]') - - ax[1].plot( - xyo_det[det_key][:, 0], - np.degrees(xyo_det[det_key][:, 2]), 'k.' - ) - ax[1].plot(xyo_i[det_key][:, 0], np.degrees(xyo_i[det_key][:, 2]), 'rx') - ax[1].grid(True) - ax[1].set_xlim(-0.5*panel.col_dim, 0.5*panel.col_dim) - ax[1].set_ylim(ome_period[0], ome_period[1]) - ax[1].set_xlabel('detector X [mm]') - ax[1].set_ylabel(r'$\omega$ [deg]') - - fig.show() - -# ============================================================================= -# %% RUN OPTIMIZATION -# ============================================================================= - -params, resd, xyo_f = calibrate_instrument_from_sx( - instr, grain_parameters, bmat, xyo_det, hkls, - ome_period=np.radians(ome_period) -) - -# define difference vectors for spot fits -for det_key, panel in instr.detectors.iteritems(): - x_diff = abs(xyo_det[det_key][:, 0] - xyo_f[det_key][:, 0]) - y_diff = abs(xyo_det[det_key][:, 1] - xyo_f[det_key][:, 1]) - ome_diff = np.degrees( - xfcapi.angularDifference( - xyo_det[det_key][:, 2], - xyo_f[det_key][:, 2]) - ) - - # filter out reflections with centroids more than - # a pixel and delta omega away from predicted value - idx_1 = np.logical_and( - x_diff <= n_pixels_tol*panel.pixel_size_col, - np.logical_and( - y_diff <= n_pixels_tol*panel.pixel_size_row, - ome_diff <= ome_tol - ) - ) - - print("INFO: Will keep %d of %d input reflections on panel %s for re-fit" - % (sum(idx_1), sum(idx_0[det_key]), det_key) - ) - - idx_new = np.zeros_like(idx_0[det_key], dtype=bool) - idx_new[np.where(idx_0[det_key])[0][idx_1]] = True - idx_0[det_key] = idx_new - -# ============================================================================= -# %% Look ok? Then proceed with refit -# ============================================================================= -# -# define difference vectors for spot fits -# for det_key, panel in instr.detectors.iteritems(): -# hkls_refit = hkls[det_key][idx_new[det_key], :] -# xyo_det_refit = xyo_det[det_key][idx_0[det_key], :] - -# update calibration crystal params -grain_parameters[:3] = params[5:8] -grain_parameters[3:6] = params[8:11] - -# reparse data -hkls_refit, xyo_det_refit, idx_0 = parse_reflection_tables( - cfg, instr, grain_id, refit_idx=idx_0 -) - -# perform refit -params2, resd2, xyo_f2 = calibrate_instrument_from_sx( - instr, grain_parameters, bmat, xyo_det_refit, hkls_refit, - ome_period=np.radians(ome_period) -) - - -# ============================================================================= -# %% plot results -# ============================================================================= - -for det_key, panel in instr.detectors.iteritems(): - fig, ax = plt.subplots(1, 2, sharex=True, sharey=False, figsize=(9, 5)) - fig.suptitle("detector %s" % det_key) - ax[0].plot(xyo_det[det_key][:, 0], xyo_det[det_key][:, 1], 'k.') - ax[0].plot(xyo_i[det_key][:, 0], xyo_i[det_key][:, 1], 'rx') - ax[0].plot(xyo_f2[det_key][:, 0], xyo_f2[det_key][:, 1], 'b+') - ax[0].grid(True) - ax[0].axis('equal') - ax[0].set_xlim(-0.5*panel.col_dim, 0.5*panel.col_dim) - ax[0].set_ylim(-0.5*panel.row_dim, 0.5*panel.row_dim) - ax[0].set_xlabel('detector X [mm]') - ax[0].set_ylabel('detector Y [mm]') - - ax[1].plot( - xyo_det[det_key][:, 0], - np.degrees(xyo_det[det_key][:, 2]), 'k.' - ) - ax[1].plot(xyo_i[det_key][:, 0], np.degrees(xyo_i[det_key][:, 2]), 'rx') - ax[1].plot(xyo_f2[det_key][:, 0], np.degrees(xyo_f2[det_key][:, 2]), 'b+') - ax[1].grid(True) - ax[1].set_xlim(-0.5*panel.col_dim, 0.5*panel.col_dim) - ax[1].set_ylim(ome_period[0], ome_period[1]) - ax[1].set_xlabel('detector X [mm]') - ax[1].set_ylabel(r'$\omega$ [deg]') - - ax[0].axis('tight') - - fig.show() - - -# ============================================================================= -# %% output results -# ============================================================================= - -# update calibration crystal params -grain_parameters[:3] = params2[5:8] -grain_parameters[3:6] = params2[8:11] - -calibration_dict = { - 'grain_id': grain_id, - 'inv_stretch': grain_parameters[6:].tolist(), - 'orientation': grain_parameters[:3].tolist(), - 'position': grain_parameters[3:6].tolist(), -} - -# write out -output_fname = cfg.instrument.parameters -if not overwrite_cfg: - output_fname = 'new-instrument.yml' -instr.write_config(filename=output_fname, - calibration_dict=calibration_dict) diff --git a/scripts/calibration_mockup.py b/scripts/calibration_mockup.py deleted file mode 100644 index a5065c3ed..000000000 --- a/scripts/calibration_mockup.py +++ /dev/null @@ -1,469 +0,0 @@ -import os - -import glob - -import matplotlib.gridspec as gridspec -from matplotlib import pyplot as plt - -import numpy as np - -from scipy.optimize import leastsq, least_squares -try: - import dill as cpl -except(ImportError): - import pickle as cpl - -from skimage.exposure import equalize_adapthist, rescale_intensity -from skimage.morphology import disk, binary_erosion - -import yaml - -from hexrd import imageseries -from hexrd import instrument -from hexrd import material -from hexrd.matrixutil import findDuplicateVectors -from hexrd.fitting import fitpeak -from hexrd.rotations import \ - angleAxisOfRotMat, \ - angles_from_rmat_xyz, make_rmat_euler, \ - RotMatEuler -from hexrd.transforms.xfcapi import mapAngle - - -# plane data -def load_pdata(cpkl, key): - with open(cpkl, "rb") as matf: - mat_list = cpl.load(matf) - return dict(list(zip([i.name for i in mat_list], mat_list)))[key].planeData - - -# images -def load_images(yml): - return imageseries.open(yml, format="image-files") - - -# instrument -def load_instrument(yml): - with open(yml, 'r') as f: - icfg = yaml.safe_load(f) - return instrument.HEDMInstrument(instrument_config=icfg) - - -# build a material on the fly -def make_matl(mat_name, sgnum, lparms, hkl_ssq_max=50): - 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 - - -# %% -# ============================================================================= -# CLASSES -# ============================================================================= - - -class InstrumentCalibrator(object): - def __init__(self, *args): - assert len(args) > 0, \ - "must have at least one calibrator" - self._calibrators = args - self._instr = self._calibrators[0].instr - - @property - def instr(self): - return self._instr - - @property - def calibrators(self): - return self._calibrators - - # ========================================================================= - # METHODS - # ========================================================================= - - def run_calibration(self, use_robust_optimization=False): - """ - FIXME: only coding serial powder case to get things going. Will - eventually figure out how to loop over multiple calibrator classes. - All will have a reference the same instrument, but some -- like single - crystal -- will have to add parameters as well as contribute to the RHS - """ - calib_class = self.calibrators[0] - - obj_func = calib_class.residual - - data_dict = calib_class._extract_powder_lines() - - # grab reduced optimizaion parameter set - x0 = self._instr.calibration_parameters[ - self._instr.calibration_flags - ] - - resd0 = obj_func(x0, data_dict) - - if use_robust_optimization: - oresult = least_squares( - obj_func, x0, args=(data_dict, ), - method='trf', loss='soft_l1' - ) - x1 = oresult['x'] - else: - x1, cox_x, infodict, mesg, ierr = leastsq( - obj_func, x0, args=(data_dict, ), - full_output=True - ) - resd1 = obj_func(x1, data_dict) - - delta_r = sum(resd0**2)/float(len(resd0)) - \ - sum(resd1**2)/float(len(resd1)) - - if delta_r > 0: - print(('OPTIMIZATION SUCCESSFUL\nfinal ssr: %f' % sum(resd1**2))) - print(('delta_r: %f' % delta_r)) - # self.instr.write_config(instrument_filename) - else: - print('no improvement in residual!!!') - - -# %% -class PowderCalibrator(object): - def __init__(self, instr, plane_data, img_dict, - tth_tol=None, eta_tol=0.25, - pktype='pvoigt'): - assert list(instr.detectors.keys()) == list(img_dict.keys()), \ - "instrument and image dict must have the same keys" - self._instr = instr - self._plane_data = plane_data - self._img_dict = img_dict - - # for polar interpolation - self._tth_tol = tth_tol or np.degrees(plane_data.tThWidth) - self._eta_tol = eta_tol - - # for peak fitting - # ??? fitting only, or do alternative peak detection? - self._pktype = pktype - - @property - def instr(self): - return self._instr - - @property - def plane_data(self): - return self._plane_data - - @property - def img_dict(self): - return self._img_dict - - @property - def tth_tol(self): - return self._tth_tol - - @tth_tol.setter - def tth_tol(self, x): - assert np.isscalar(x), "tth_tol must be a scalar value" - self._tth_tol = x - - @property - def eta_tol(self): - return self._eta_tol - - @eta_tol.setter - def eta_tol(self, x): - assert np.isscalar(x), "eta_tol must be a scalar value" - self._eta_tol = x - - @property - def pktype(self): - return self._pktype - - @pktype.setter - def pktype(self, x): - """ - currently only 'pvoigt' or 'gaussian' - """ - assert isinstance(x, str), "tth_tol must be a scalar value" - self._pktype = x - - def _interpolate_images(self): - """ - returns the iterpolated powder line data from the images in img_dict - - ??? interpolation necessary? - """ - return self.instr.extract_line_positions( - self.plane_data, self.img_dict, - tth_tol=self.tth_tol, eta_tol=self.eta_tol, - npdiv=2, collapse_eta=False, collapse_tth=False, - do_interpolation=True) - - def _extract_powder_lines(self): - """ - return the RHS for the instrument DOF and image dict - - The format is a dict over detectors, each containing - - [index over ring sets] - [index over azimuthal patch] - [xy_meas, tth_meas, tth_ref, eta_ref] - - FIXME: can not yet handle tth ranges with multiple peaks! - """ - # ideal tth - tth_ideal = self.plane_data.getTTh() - tth0 = [] - for idx in self.plane_data.getMergedRanges()[0]: - if len(idx) > 1: - eqv, uidx = findDuplicateVectors(np.atleast_2d(tth_ideal[idx])) - if len(uidx) > 1: - raise NotImplementedError("can not handle multipeak yet") - else: - # if here, only degenerate ring case - uidx = idx[0] - else: - uidx = idx[0] - tth0.append(tth_ideal[uidx]) - - powder_lines = self._interpolate_images() - - # GRAND LOOP OVER PATCHES - rhs = dict.fromkeys(self.instr.detectors) - for det_key, panel in self.instr.detectors.items(): - rhs[det_key] = [] - for i_ring, ringset in enumerate(powder_lines[det_key]): - tmp = [] - for angs, intensities in ringset: - tth_centers = np.average( - np.vstack([angs[0][:-1], angs[0][1:]]), - axis=0) - eta_ref = angs[1] - int1d = np.sum(np.array(intensities).squeeze(), axis=0) - """ - DARREN: FIT [tth_centers, intensities[0]] HERE - - RETURN TTH0 - rhs.append([tth0, eta_ref]) - """ - p0 = fitpeak.estimate_pk_parms_1d( - tth_centers, int1d, self.pktype - ) - - p = fitpeak.fit_pk_parms_1d( - p0, tth_centers, int1d, self.pktype - ) - # !!! this is where we can kick out bunk fits - tth_meas = p[1] - center_err = abs(tth_meas - tth0[i_ring]) - if p[0] < 0.1 or center_err > np.radians(self.tth_tol): - continue - xy_meas = panel.angles_to_cart([[tth_meas, eta_ref], ]) - - # FIXME: distortion kludge - if panel.distortion is not None: - xy_meas = panel.distortion[0]( - xy_meas, - panel.distortion[1], - invert=True - ) - # cat results - tmp.append( - np.hstack( - [xy_meas.squeeze(), - tth_meas, - tth0[i_ring], - eta_ref] - ) - ) - rhs[det_key].append(np.vstack(tmp)) - rhs[det_key] = np.vstack(rhs[det_key]) - return rhs - - def residual(self, reduced_params, data_dict): - """ - """ - - # first update instrument from input parameters - full_params = self.instr.calibration_parameters - full_params[self.instr.calibration_flags] = reduced_params - self.instr.update_from_parameter_list(full_params) - - # build residual - resd = [] - for det_key, panel in self.instr.detectors.items(): - pdata = np.vstack(data_dict[det_key]) - if len(pdata) > 0: - calc_xy = panel.angles_to_cart(pdata[:, -2:]) - - # FIXME: distortion kludge - if panel.distortion is not None: - calc_xy = panel.distortion[0]( - calc_xy, - panel.distortion[1], - invert=True - ) - - resd.append( - (pdata[:, :2].flatten() - calc_xy.flatten()) - ) - else: - continue - return np.hstack(resd) - - -# %% -# ============================================================================= -# USER INPUT -# ============================================================================= - -# dirs -working_dir = '/Users/Shared/APS/PUP_AFRL_Feb19' -image_dir = os.path.join(working_dir, 'image_data') - -samp_name = 'ceria_cal' -scan_number = 0 - -raw_data_dir_template = os.path.join( - image_dir, - 'raw_images_%s_%06d-%s.yml' -) - -instrument_filename = 'Hydra_Feb19.yml' - -# !!! make material on the fly as pickles are broken -matl = make_matl('ceria', 225, [5.411102, ]) - -# tolerances for patches -tth_tol = 0.2 -eta_tol = 2. - -tth_max = 11. - -# peak fit type -# pktype = 'gaussian' -pktype = 'pvoigt' - -# Try it -use_robust_optimization = False - - -# %% -# ============================================================================= -# IMAGESERIES -# ============================================================================= - -# for applying processing options (flips, dark subtretc...) -PIS = imageseries.process.ProcessedImageSeries - -# load instrument -#instr = load_instrument(os.path.join(working_dir, instrument_filename)) -instr = load_instrument(instrument_filename) -det_keys = list(instr.detectors.keys()) - -# hijack panel buffer -edge_buff = 5 -for k, v in list(instr.detectors.items()): - try: - mask = np.load(mask_template % k) - except(NameError): - mask = np.ones((v.rows, v.cols), dtype=bool) - mask[:edge_buff, :] = False - mask[:, :edge_buff] = False - mask[-edge_buff:, :] = False - mask[:, -edge_buff:] = False - v.panel_buffer = mask - - -# grab imageseries filenames -file_names = glob.glob( - raw_data_dir_template % (samp_name, scan_number, '*') -) -check_files_exist = [os.path.exists(file_name) for file_name in file_names] -if not np.all(check_files_exist): - raise RuntimeError("files don't exist!") - -img_dict = dict.fromkeys(det_keys) -for k, v in img_dict.items(): - file_name = glob.glob( - raw_data_dir_template % (samp_name, scan_number, k) - ) - ims = load_images(file_name[0]) - img_dict[k] = ims[0] - -# plane data munging -matl.beamEnergy = instr.beam_energy -plane_data = matl.planeData -if tth_tol is not None: - plane_data.tThWidth = np.radians(tth_tol) - -if tth_max is not None: - plane_data.exclusions = None - plane_data.tThMax = np.radians(tth_max) - - -# %% -# ============================================================================= -# SETUP CALIBRATION -# ============================================================================= - -# first, add tilt calibration -# !!! This needs to come from the GUI input somehow -rme = RotMatEuler(np.zeros(3), 'xyz', extrinsic=True) -instr.tilt_calibration_mapping = rme - -# !!! this comes from GUI checkboxes -# set flags for first 2 tilt angles, translation for each panel (default) -# first three distortion params -cf = instr.calibration_flags -ii = 7 -for i in range(instr.num_panels): - cf[ii + 2] = False - cf[ii + 6:ii + 9] = True - ii += 12 -instr.calibration_flags = cf - -# powder calibrator -pc = PowderCalibrator(instr, plane_data, img_dict, - tth_tol=tth_tol, eta_tol=eta_tol) - -# make instrument calibrator -ic = InstrumentCalibrator(pc) - -ic.run_calibration() - -# %% sample plot to check fit line poistions ahead of fitting -fig, ax = plt.subplots(2, 2) -fig_row, fig_col = np.unravel_index(np.arange(instr.num_panels), (2, 2)) - -data_dict = pc._extract_powder_lines() - -ifig = 0 -for det_key, panel in instr.detectors.items(): - all_pts = np.vstack(data_dict[det_key]) - ''' - pimg = equalize_adapthist( - rescale_intensity(img_dict[det_key], out_range=(0, 1)), - 10, clip_limit=0.2) - ''' - pimg = np.array(img_dict[det_key], dtype=float) - pimg[~panel.panel_buffer] = np.nan - ax[fig_row[ifig], fig_col[ifig]].imshow( - pimg, - vmin=np.percentile(img_dict[det_key], 5), - vmax=np.percentile(img_dict[det_key], 90), - cmap=plt.cm.bone_r - ) - ideal_angs, ideal_xys = panel.make_powder_rings(plane_data, - delta_eta=eta_tol) - rijs = panel.cartToPixel(np.vstack(ideal_xys)) - ax[fig_row[ifig], fig_col[ifig]].plot(rijs[:, 1], rijs[:, 0], 'cx') - ax[fig_row[ifig], fig_col[ifig]].set_title(det_key) - rijs = panel.cartToPixel(all_pts[:, :2]) - ax[fig_row[ifig], fig_col[ifig]].plot(rijs[:, 1], rijs[:, 0], 'm+') - ax[fig_row[ifig], fig_col[ifig]].set_title(det_key) - ifig += 1 diff --git a/scripts/convert_instrument_config.py b/scripts/convert_instrument_config.py deleted file mode 100644 index 2cea77ccd..000000000 --- a/scripts/convert_instrument_config.py +++ /dev/null @@ -1,161 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -This module provides a cmd-line callable utility for converting older-style -instrument YAML configuration files to hexrd >=0.6 - -Examples --------- - $ python convert_instrument_config.py --help - $ python convert_instrument_config.py -o new_param.yml - -Notes ------ - - Note that the format changes include the following: - * All translation vector components now have the key `translation` rather - than `t_vec_*` for clarity. - * The `tilt_angles` for each detector are now under `tilt`, and have been - changed from extrinsic XYZ Euler angles (in radians) to the associated - rotation matrix invariants; specifically angle*axis (dimensionless). - This was done to more easily facilitate the use of different user- - specified Euler angle conventions in external interfaces. - -Todo: - * Add to hexrd/scripts -""" - -import argparse -import yaml -from hexrd.rotations import make_rmat_euler, angleAxisOfRotMat - -# ============================================================================= -# PARAMETERS -# ============================================================================= - -# top level keys -ROOT_KEYS = ['beam', 'oscillation_stage', 'detectors'] -OLD_OSCILL_KEYS = ['t_vec_s', 'chi'] -NEW_OSCILL_KEYS = ['translation', 'chi'] -OLD_TRANSFORM_KEYS = ['t_vec_d', 'tilt_angles'] -NEW_TRANSFORM_KEYS = ['translation', 'tilt'] - - -# ============================================================================= -# FUNCTIONS -# ============================================================================= - -def convert_instrument_config(old_cfg, output=None): - """ - Convert v0.5.x style instrument YAML config file to v0.6.x style. - - Parameters - ---------- - old_cfg : str - The filename of the v0.5.x style instrument YAML config. - output : str, optional - The filename to write the updated v0.6.x instrument YAML config to. - The default is None, in which case the new YAML format is printed to - stdout. - - Raises - ------ - RuntimeError - Raises exception if input config has any unrecognized keys. - - Returns - ------- - None. - - """ - icfg = yaml.safe_load(open(old_cfg, 'r')) - - new_cfg = dict.fromkeys(icfg) - - # %% first beam - new_cfg['beam'] = icfg['beam'] - - # %% next, calibration crystal if applicable - calib_key = 'calibration_crystal' - if calib_key in icfg.keys(): - new_cfg[calib_key] = icfg[calib_key] - - # %% sample stage - old_dict = icfg['oscillation_stage'] - tmp_dict = dict.fromkeys(NEW_OSCILL_KEYS) - for tk in zip(OLD_OSCILL_KEYS, NEW_OSCILL_KEYS): - tmp_dict[tk[1]] = old_dict[tk[0]] - new_cfg['oscillation_stage'] = tmp_dict - - # %% detectors - new_cfg['detectors'] = dict.fromkeys(icfg['detectors']) - det_block_keys = ['pixels', 'saturation_level', 'transform', 'distortion'] - for det_id, source_params in icfg['detectors'].items(): - new_dict = {} - for key in det_block_keys: - if key != 'transform': - try: - new_dict[key] = source_params[key] - except(KeyError): - if key == 'distortion': - continue - elif key == 'saturation_level': - new_dict[key] = 2**16 - else: - raise RuntimeError( - "unrecognized parameter key '%s'" % key - ) - else: - old_dict = source_params[key] - tmp_dict = dict.fromkeys(NEW_TRANSFORM_KEYS) - for tk in zip(OLD_TRANSFORM_KEYS, NEW_TRANSFORM_KEYS): - if tk[0] == 't_vec_d': - tmp_dict[tk[1]] = old_dict[tk[0]] - elif tk[0] == 'tilt_angles': - xyz_angles = old_dict[tk[0]] - rmat = make_rmat_euler(xyz_angles, - 'xyz', - extrinsic=True) - phi, n = angleAxisOfRotMat(rmat) - tmp_dict[tk[1]] = (phi*n.flatten()).tolist() - new_dict[key] = tmp_dict - new_cfg['detectors'][det_id] = new_dict - - # %% dump new file - if output is None: - print(new_cfg) - else: - yaml.dump(new_cfg, open(output, 'w')) - return - - -# ============================================================================= -# EXECUTION -# ============================================================================= - -if __name__ == '__main__': - parser = argparse.ArgumentParser( - description="convert a v0.5.x instrument config to v0.6.x format" - ) - - parser.add_argument( - 'instrument_cfg', - help="v0.5 style instrument config YAML file" - ) - - parser.add_argument( - '-o', '--output-file', - help="output file name", - type=str, - default="" - ) - - args = parser.parse_args() - - old_cfg = args.instrument_cfg - output_file = args.output_file - - if len(output_file) == 0: - output_file = None - - convert_instrument_config(old_cfg, output=output_file) diff --git a/scripts/pinhole_2theta.py b/scripts/pinhole_2theta.py deleted file mode 100644 index 6a8f84ab7..000000000 --- a/scripts/pinhole_2theta.py +++ /dev/null @@ -1,836 +0,0 @@ -""" -Calculation of x-ray diffraction scattering angle for PXRDIP and TARDIS. - -Calculates the nominal, sample-weighted, and pinhole-weighted twotheta - scattering angle for PXRDIP and TARDIS diffraction experiments. - -Reference: -[1] J. R. Rygg et al, "X-ray diffraction at the National Ignition Facility," - Rev. Sci. Instrum. 91, 043902 (2020). https://doi.org/10.1063/1.5129698 - -file: pinhole_2theta.py -revisions: - 2022-02-02, v1.5: refactor for standalone script, additional comments, - prepare for commit to https://github.com/HEXRD/hexrd - 2021-11-18, v1.4: option to increase twotheta precision for output table - 2020-12-21, v1.3: additional metadata in report file and filename - 2020-09-21, v1.2: added ReportFileWriter to automate file writing - 2020-08-22, v1.1: minor updates - 2020-08-09, v1.0: migrated from earlier scripts -""" -# written in python 3.7; not tested with earlier versions. -__author__ = "J. R. Rygg" -__version__ = "1.5" - -# header of text file report -HEADER = ( -'This file contains tables of the evaluated pinhole twotheta correction values\n' -' described in section III.C.2 of ref (1) for parameters defined below.\n' -'\n' -' (1) J. R. Rygg et al, "X-ray Diffraction at the National Ignition Facility,"\n' -' Rev. Sci. Instrum. 91, 043902 (2020). https://doi.org/10.1063/1.5129698\n' -'\n[definitions]\n' -'twotheta_n: nominal twotheta, Ref (1) Eq. 5.\n' -'twotheta_p: pinhole twotheta, Ref (1) section III.C.2 and Eq 6.\n' -'twotheta_JHE: J. Eggert approximation of twotheta_p, "In my PXRDIP and\n' -' TARDIS Igor analysis package, I approximated this correction by\n' -' assuming that the pinhole diffraction arises from the center of the\n' -' pinhole wall opposite the image-plate pixel."\n' -) - -# imports -import datetime -import io -import os.path as osp - -import numpy as np -from numpy import pi, sin, cos, tan, radians - -import matplotlib as mpl -import matplotlib.pyplot as plt - -DPI_SCREEN = 150 - -# reset mpl rc to defaults and adjust settings for image generation -mpl.rcdefaults() -mpl.rcParams.update({'backend': 'QtAgg' }) -mpl.rc('font',**{'family':'sans-serif', 'sans-serif':['Arial'], 'size':12}) -mpl.rc('mathtext', default='regular') # same font for mathtext as normal text -mpl.rc('image', cmap='viridis', origin='lower') # -mpl.rc('lines', linewidth=0.75) # -mpl.rc('xtick', top=True, bottom=True) -mpl.rc('ytick', left=True, right=True) - -# linear attenuation coefficient [um^-1] for pinhole, xsf pairs -MU_P_LOOKUP = { # (pinhole_mat, xs_mat): linear attenuation coefficient [um^-1] - ('Ta', 'Fe'): 1/2.33, ('Ta', 'Cu'): 1/4.11, ('Ta', 'Ge'): 1/2.69, # 1/mic - ('W', 'Fe'): 1/1.94, ('W', 'Cu'): 1/3.42, ('W', 'Ge'): 1/2.25, # 1/mic - ('Pt', 'Fe'): 1/1.50, ('Pt', 'Cu'): 1/2.63, ('Pt', 'Ge'): 1/4.38, # 1/mic -} - -# baseline parameters -OMEGA_BASELINE=dict(mat_x='Cu', r_x=24.24, alpha=45, phi_x=0, - mat_p='Ta', d_p=0.3, h_p=0.075, centered_ph=True) -EP_BASELINE =dict(mat_x='Cu', r_x=24.24, alpha=22.5, phi_x=0, - mat_p='Ta', d_p=0.3, h_p=0.075, centered_ph=True) -NIF_BASELINE =dict(mat_x='Ge', r_x=32, alpha=30, phi_x=0, - mat_p='Ta', d_p=0.4, h_p=0.1, centered_ph=True) - -# output path for image and report text file -OUTPUT_PATH = osp.dirname(__file__) - -# ============================ Helper Functions =============================== -def zenith(vv, v0=(0,0,1), uom='radians'): - """Return zenith angle between vector array vv and reference vector(s) v0. - - Sometimes known as the polar angle, inclination angle, or colatitude. - expects the 3-vector to be contained along last axis. - """ - try: - np.broadcast(vv, v0) - except ValueError: # if vv and v0 as given cannot be broadcast to eachother - shapev, shape0 = np.shape(vv), np.shape(v0) - shapev_new = shapev[:-1] + (1,) * (len(shape0) - 1) + (shapev[-1],) - shape0_new = (1,) * (len(shapev) - 1) + shape0 - vv = np.reshape(vv, shapev_new) - v0 = np.reshape(v0, shape0_new) - - # normalization factor along last axis - norm = np.linalg.norm(vv, axis=-1) * np.linalg.norm(v0, axis=-1) - - # np.clip eliminates rounding errors for (anti)parallel vectors - out = np.arccos(np.clip(np.sum(vv * v0, -1) / norm, -1, 1)) - if 'deg' in uom: - return out * 180 / pi - return out - - -def azimuth(vv, v0=None, v1=None, uom='radians'): - """Return azimuthal angle btwn vv and v0, with v1 defining phi=0.""" - if v0 is None: v0 = np.array((0, 0, 1), dtype=float) # z unit vector - if v1 is None: v1 = np.array((1, 0, 0), dtype=float) # x unit vector - - with np.errstate(divide='ignore', invalid='ignore'): - n0 = np.cross(v0, v1) - n0 /= np.linalg.norm(n0, axis=-1)[...,np.newaxis] - nn = np.cross(v0, vv) - nn /= np.linalg.norm(nn, axis=-1)[...,np.newaxis] - - azi = np.arccos(np.sum(nn * n0, -1)) - if len(np.shape(azi)) > 0: - azi[np.dot(vv, n0) < 0] *= -1 - azi[np.isnan(azi)] = 0 # arbitrary angle where vv is (anti)parallel to v0 - else: - if np.isnan(azi): - return 0 - elif np.dot(vv, v0) < 1 and azi > 0: - azi *= -1 - - if 'deg' in uom: - return azi * 180 / pi - return azi - - -def calc_twotheta_nominal(alpha, beta, phi_x=0, phi_d=0): # Eq. (5) - """Nominal twotheta [radians] given alpha, beta, phi_x, phi_d. - - optional: alpha, phi_x taken from first argument if is XraySource instance - """ - if isinstance(alpha, XraySource): - alpha, phi_x = alpha.alpha, alpha.phi - - # Rygg2020, Eq. (5) - cos_2qn = cos(alpha)*cos(beta) - sin(alpha)*sin(beta) * cos(phi_d - phi_x) - return np.arccos(np.clip(cos_2qn, -1, 1)) - - -def calc_twotheta_sample_minus_nominal(qq_n, beta, r_x, z_s): # Eq. (8) - """Return (sample twotheta - nominal twotheta) [radians]. - - See Rygg2020, Eq. (8) - """ - return np.arctan(sin(qq_n) / (r_x/z_s * cos(beta) - cos(qq_n))) - - -def calc_twotheta_JHE(X, P, D): - """J. Eggert approximation to pinhole twotheta. - - Quote from J. Eggert describing approximation: "In my PXRDIP and TARDIS - Igor analysis package, I approximated this correction by assuming that - the pinhole diffraction arises from the center of the pinhole wall - opposite the image-plate pixel." - """ - v_x = X.v_x if isinstance(X, XraySource) else X - d_p = P.diam if isinstance(P, Pinhole) else P - v_d = D.v_d if isinstance(D, Detector) else D - - b = d_p / 2 # impact parameter - dw = np.zeros_like(v_d) # offset vector of sample interaction region - dw[...,0] = b * cos(D.phi + pi) - dw[...,1] = b * sin(D.phi + pi) - - return zenith(v_d - dw, v0=(v_x + dw)) - - -def calc_twotheta_pinhole(X, P, D, Np=120, showdetailview=False): - """Return pinhole twotheta [rad] and effective scattering volume [mm3]. - - X: XraySource instance, or position tuple - P: Pinhole instance, or (material, diameter, thickness) tuple - D: Detector insteance, or tuple of (r_d, beta, phi_d) - Np: number of pinhole phi elements for integration - showdetailview: show detailed report, and plots if one detector element - """ - #------ Determine geometric parameters ------ - if isinstance(X, XraySource): - r_x, alpha, phi_x = X.r, X.alpha, X.phi - mat_x = X.mat - else: # assume tuple of r_x, alpha, phi_x - r_x, alpha, phi_x = X - mat_x = "Cu" - if isinstance(P, Pinhole): - mat_p, d_p, h_p = P.mat, P.d, P.h - else: # assume tuple of (material, diameter, thickness) - mat_p, d_p, h_p = P - - mu_p = MU_P_LOOKUP.get((mat_p, mat_x), 1/3.0) # attenuation coefficient - mu_p = 1000 * mu_p # convert to mm**-1 - - if isinstance(D, Detector): - r_d, beta, phi_d = D.r_d, D.beta, D.phi_d - elif len(D) == 3: # assume tuple of r_d, beta, phi_d - r_d, beta, phi_d = D - else: # handle other inputs - print("unexpected form for 'D' in calc_twotheta_pinhole") - - # reshape so D grids are on axes indices 2 and 3 [1 x 1 x Nu x Nv] - r_d = np.atleast_2d(r_d)[None, None, :, :] - beta = np.atleast_2d(beta)[None, None, :, :] - phi_d = np.atleast_2d(phi_d)[None, None, :, :] - - #------ define pinhole grid ------ - Nz = max(3, int(Np * h_p / (pi * d_p))) # approximately square pinhole surface elements - dphi = 2 * pi / Np # [rad] phi interval - dl = d_p * dphi # [mm] azimuthal distance increment - dz = h_p / Nz # [mm] axial distance increment - dA = dz * dl # [mm^2] area element - dV_s = dA * mu_p**-1 # [mm^3] volume of surface element - dV_e = dl * mu_p**-2 # [mm^3] volume of edge element - - phi_vec = np.arange(dphi/2, 2*pi, dphi) - z_vec = np.arange(-h_p/2 - dz/2, h_p/2 + dz/1.999, dz) # includes elements for X and D edges - z_vec[0] = -h_p/2 # X-side edge (negative z) - z_vec[-1] = h_p/2 # D-side edge (positive z) - phi_i, z_i = np.meshgrid(phi_vec, z_vec) # [Nz x Np] - phi_i = phi_i[:, :, None, None] # [Nz x Np x 1 x 1] - z_i = z_i[:, :, None, None] # axes 0,1 => P; axes 2,3 => D - - #------ calculate twotheta_i [a.k.a. qq_i], for each grid element ------ - bx, bd = (d_p / (2 * r_x), d_p / (2 * r_d)) - sin_a, cos_a, tan_a = sin(alpha), cos(alpha), tan(alpha) - sin_b, cos_b, tan_b = sin(beta), cos(beta), tan(beta) - sin_phii, cos_phii = sin(phi_i), cos(phi_i) - cos_dphi_x = cos(phi_i - phi_x + pi) # [Nz x Np x Nu x Nv] - cos_dphi_d = cos(phi_i - phi_d + pi) - - arctan2 = np.arctan2 - alpha_i = arctan2(np.sqrt(sin_a**2 + 2*bx*sin_a*cos_dphi_x + bx**2), cos_a + z_i/r_x) - beta_i = arctan2(np.sqrt(sin_b**2 + 2*bd*sin_b*cos_dphi_d + bd**2), cos_b - z_i/r_d) - phi_xi = arctan2(sin_a*sin(phi_x) - bx*sin_phii, sin_a*cos(phi_x) - bx*cos_phii) - phi_di = arctan2(sin_b*sin(phi_d) - bd*sin_phii, sin_b*cos(phi_d) - bd*cos_phii) - - arg = cos(alpha_i) * cos(beta_i) - sin(alpha_i) * sin(beta_i) * cos(phi_di - phi_xi) - qq_i = np.arccos(np.clip(arg, -1, 1)) # scattering angle for each P to each D - - #------ calculate effective volumes: 1 (surface), 2 (Xedge), 3 (Dedge) ----- - sec_psi_x = 1 / (sin_a * cos_dphi_x) - sec_psi_d = 1 / (sin_b * cos_dphi_d) - sec_alpha = 1 / cos_a - sec_beta = 1 / cos_b - tan_eta_x = np.where(cos_dphi_x[0] <= 0, 0, cos_a * cos_dphi_x[0]) - tan_eta_d = np.where(cos_dphi_d[-1] <= 0, 0, cos_b * cos_dphi_d[-1]) - - V_i = dV_s / (sec_psi_x + sec_psi_d) # [mm^3] - V_i[0] = dV_e / (sec_psi_d[0] * (sec_alpha + sec_psi_d[0] * tan_eta_x)) # X-side edge (z = -h_p / 2) - V_i[-1] = dV_e / (sec_psi_x[-1] * (sec_beta + sec_psi_x[-1] * tan_eta_d)) # D-side edge (z = +h_p / 2) - - #------ visibility of each grid element ------ - is_seen = np.logical_and(z_i > h_p/2 - d_p/tan_b * cos_dphi_d, # pinhole surface - z_i < -h_p/2 + d_p/tan_a * cos_dphi_x) - is_seen[0] = np.where(h_p/d_p * tan_b < cos_dphi_d[0], 1, 0) # X-side edge - is_seen[-1] = np.where(h_p/d_p * tan_a < cos_dphi_x[-1], 1, 0) # D-side edge - - #------ weighted sum over elements to obtain average ------ - V_i *= is_seen # zero weight to elements with no view of both X and D - V_p = np.nansum(V_i, axis=(0,1)) # [Nu x Nv] <= detector - qq_p = np.nansum(V_i * qq_i, axis=(0,1)) / V_p # [Nu x Nv] <= detector - - if not showdetailview: - return qq_p, V_p - - #====== Print and plot detailed view ====== - print("=== Pinhole.calc_twotheta_pinhole: detailed view ===") - print("{} PH with (d={:.3f} mm, h={:.3f} mm, mu_p=1/({:.2f} um))".format(mat_p, d_p, h_p, 1000/mu_p)) - print("{} XS at (r={:.1f} mm, alpha={:.1f} deg, phi={:.0f})".format(mat_x, float(r_x), np.degrees(float(alpha)), np.degrees(float(phi_x)))) - print("number of detector elements =", np.size(beta)) - print("number of pinhole elements =", np.size(phi_i)) - print("number of scattering calculations =", np.size(qq_i)) - if np.size(beta) == 1: - print("detector pixel at (r={:.1f} mm, beta={:.0f} deg, phi={:.0f})".format(float(r_d), np.degrees(float(beta)), np.degrees(float(phi_d)))) - qq_n = np.squeeze(calc_twotheta_nominal(alpha, beta, phi_x, phi_d)) - qq_p = np.squeeze(qq_p) -# qq_JHE = np.squeeze(calc_twotheta_JHE(X, d_p, D)) - - print("nominal twotheta = {:7.4f} deg".format(np.degrees(qq_n))) - print("pinhole twotheta = {:7.4f} deg".format(np.degrees(qq_p))) - print("twotheta correc. = {:7.4f} deg".format(np.degrees(qq_p - qq_n))) - print("max 2q deviation = {:7.4f} deg".format(np.degrees(np.max(np.squeeze(qq_i)) - qq_n))) - -# qq_JHE = self.calc_twotheta_pinhole_JHE(X, (r_d, beta, phi_d)) -# print("JHE approx = {:7.4f} deg".format(np.degrees(qq_JHE - qq_n))) - - print("signal volume (total) = {:.2g} um3".format(1e9*np.sum(V_i))) - print("signal volume (total) = {:.1f}% of surface ref volume (pi d h / mu)".format(100*np.sum(V_i)/(pi*d_p*h_p/mu_p))) - print("signal contribution from surface = {:.1f}%".format(100*np.sum(V_i[1:-1,:])/np.sum(V_i))) - - extent = (0,360,-h_p/2,h_p/2) - fig = plt.figure() - - ax1 = plt.subplot(311) - im = ax1.imshow(np.degrees(np.squeeze(qq_i[1:-1,:] - qq_n)), aspect='auto', - vmin=-0.6, vmax=0.6, cmap='coolwarm',extent=extent) - plt.colorbar(im, label=r"$2\theta_p - 2\theta_n$ [deg]") - ax1.plot(180/pi*(np.squeeze(phi_d)+pi), 0, 'wD') -# cb.set_label("2theta_p - 2theta_n [deg]") - - ax2 = plt.subplot(312, sharex=ax1, sharey=ax1) - im = ax2.imshow(1e9*np.squeeze(V_i), aspect='auto',extent=extent) -# im = ax.imshow(np.squeeze(V_i[1:-1,:]), aspect='auto',extent=extent) - plt.colorbar(im, label=r"$V_i$ [um$^3$]") - ax2.plot(180/pi*(np.squeeze(phi_d)+pi), 0, 'wD') - - ax3 = plt.subplot(313, sharex=ax1, sharey=ax1, xlabel="phi_p") - im = ax3.imshow(np.squeeze(is_seen), aspect='auto', extent=extent) -# im = ax.imshow(np.squeeze(is_seen[1:-1,:]), aspect='auto', extent=extent) -# print(np.shape(phi_vec), np.shape(h_p/2 - d_p/tan_b * cos_dphi_d)) - ax3.plot(180/pi*phi_vec, np.squeeze(h_p/2 - d_p/tan_b * cos(phi_vec - phi_d + pi))) - ax3.plot(180/pi*phi_vec, np.squeeze(-h_p/2 + d_p/tan_a * cos(phi_vec - phi_x + pi))) - ax3.plot(180/pi*(np.squeeze(phi_d)+pi), 0, 'wD') - ax3.axis(extent) - fig.tight_layout() - - print("=== /end detail view ===") - return qq_p, V_p - - -# ============================ Utility Classes =============================== -class XraySource(): - """X-ray source (XS), e.g. due to He-alpha emission of a mid-Z metal. - - Parameters - ---------- - mat : backlighter foil material (e.g. Fe, Cu, Ge) - r : distance of XS from origin (i.e. center of pinhole) [mm] - alpha : zenith angle of XS from (negative) pinhole axis - phi : azimuthal angle of XS around pinhole axis - indegrees : whether the input angles are in degrees, radians if not True - - Attributes - ---------- - alpha, phi : stored in radians. degree-values can be retrieved as alpha_deg - v_x : emission vector toward origin - """ - - def __init__(self, mat='Cu', r=24.14, alpha=45, phi=0, indegrees=True): - self.mat = mat - self.r = r - if indegrees: - self.alpha_deg, self.phi_deg = alpha, phi - self.alpha = radians(alpha) - self.phi = radians(phi) - else: - self.alpha, self.phi = alpha, phi - self.alpha_deg = np.degrees(alpha) - self.phi_deg = np.degrees(phi) - - sin_alpha, cos_alpha = sin(self.alpha), cos(self.alpha) - sin_phi = np.around(sin(self.phi), 8) # np.around to handle rounding error - cos_phi = np.around(cos(self.phi), 8) - - self.v_x = r * np.array((sin_alpha * -cos_phi, # assumes XS position at -z - sin_alpha * -sin_phi, - cos_alpha)) - - self.M = self.get_M() - - def get_M(self): - if self.phi != 0: - print("Non phi_x=0 XraySource Matrix not implemented") - return np.matrix(np.eye(3)) - - cosa, sina = cos(self.alpha), sin(self.alpha) - return np.array(((cosa, 0, -sina), # transformation matrix - (0, 1, 0), - (sina, 0, cosa))) - - def transform_to_box_coords(self, twotheta, phi=0): - """Return (beta, phi_d) for given (twotheta_n, phi).""" - input_shape = np.shape(twotheta) - twotheta = np.ravel(twotheta) - phi = np.ravel(phi) - cosq, sinq = np.cos(twotheta), np.sin(twotheta) - cosp, sinp = np.cos(phi), np.sin(phi) - - # xs direct image coords, and transformed to box coords - v_1 = np.squeeze(np.transpose(np.dstack((sinq * cosp, sinq * sinp, cosq)))) # words for 1d array of vectors - v_2 = np.transpose(np.dot(self.M, v_1)) - - beta = np.reshape(zenith(v_2), input_shape) - phi_d = np.reshape(azimuth(v_2), input_shape) - - return beta, phi_d - - -class Pinhole(): - """Cylindrical aperture in an opaque material. - - The diagnostic coordinate system (DCS) uses the pinhole center as the - origin, and the pinhole axis as the z-axis - - Parameters - ---------- - mat : pinhole substrate material (e.g. Ta, W, Pt) - diam : pinhole aperture diameter [mm] - thick : pinhole substrate thickness [mm] - pos : pinhole center position; optional, default (0,0,0) - normal : pinhole normal unit vector; optional, default (0,0,1) - - Attributes - ---------- - all __init__ parameters are saved as attributes - d, h : aliases for diam and thick, respectively - critical_angle : critical angle at which pinhole is effectively closed [rad] - """ - def __init__(self, mat='Ta', diam=0.4, thick=0.1, - pos=(0,0,0), normal=(0,0,1), **kws): - self.mat = mat - self.diam = self.d = diam - self.thick = self.h = thick - self.pos = np.array(pos) - self.normal = np.array(normal) - - @property - def critical_angle(self): - return np.arctan(self.diam / self.thick) - - -class Detector(): - """A detector single element or array of elements.""" - PXRDIP_LWH = (75, 50, 25) # [mm] length, width, and offset of pxrdip box - - def __init__(self, v_d, id=None): - """v_d contains the 3-vector(s) of the detector elements in xyz.""" - self.id = id - self.v_d = v_d - self.r_d = np.sqrt(np.sum(self.v_d**2,axis=2)) # nominal beta of each pixel - self.beta = zenith(self.v_d) # nominal beta of each pixel - self.phi_d = azimuth(self.v_d) # azimuthal angle of each pixel - - @classmethod - def pxrdip_from_id(cls, plate_id='D', centered_ph=True, shape=None): - """Generate an element array corresponding to a given PXRDIP plate.""" - plate_id = plate_id.upper()[0] - L, W, _ = cls.PXRDIP_LWH # width and length of PXRDIP image plates - if plate_id == "B": - L = W # square back plate - - if shape is None: # default to 1 element per mm - shape = (int(W) + 1, int(L) + 1) - N1, N2 = shape - - # p,q in image plate coordinate system (IPCS) - p = np.linspace(0, W, N1).reshape(N1, 1) - W / 2 - q = np.linspace(0, L, N2).reshape(1, N2) - - # convert to diagnostic coordinate system (DCS) - v_d = np.empty((N1, N2, 3)) # placeholder for detector pixel positions - if plate_id == "B": - q -= L / 2 - - v_d[..., 0] = p - v_d[..., 1] = q - v_d[..., 2] = cls.PXRDIP_LWH[0] - else: - offset_sign = (-1, 1)[int(plate_id in ('R', 'U'))] - ordering = (1, -1)[int(plate_id in ('L', 'U'))] - i_offset = int(plate_id in ('L', 'R')) - i_p = 1 - i_offset # v_d index corresponding to "p" direction - - v_d[..., i_offset] = offset_sign * W / 2 - v_d[..., i_p] = p[::ordering] - v_d[..., 2] = q - - obj = cls(v_d, plate_id) - obj.plate_id = plate_id - obj.centered_ph = centered_ph - return obj - - @classmethod - def pxrdip_from_twothetaphi(cls, xs, twotheta_n=None, phi=None, centered_ph=True): - """Generate element array for given twotheta,phi on PXRDIP box.""" - L, _, H = cls.PXRDIP_LWH # mm, length of PXRDIP box, height of pinhole vs D plate - if twotheta_n is None and phi is None: - qq = np.arange(10, xs.alpha_deg + 89, 1) # zenith every 1 degrees - pp = np.arange(0, 180 + 1, 5) # azimuth every 5 degrees (calculate half since assuming symmetric so far) - twotheta_n, phi = np.meshgrid(qq, pp) - - twotheta_n = np.radians(twotheta_n) - phi = np.radians(phi) - - beta, phi_d = xs.transform_to_box_coords(twotheta_n, phi) - sinb, cosb = sin(beta), cos(beta) - - # deduce detector element distances, r_d - rmax = np.sqrt(L**2 + 2 * H**2) - back_betamax = np.arccos(L / rmax) - down_betamin = np.arccos(L / np.sqrt(L**2 + 1 * H**2)) - r_back = np.where(beta < back_betamax, L / cosb, rmax) - r_drul = np.where(np.logical_and(down_betamin < beta, beta < np.pi/2), - H / (sinb * cos(radians((180/pi*phi_d+45)%90 - 45))), rmax) - r_d = np.minimum(r_back, r_drul) - - x = r_d * sinb * cos(phi_d) - y = r_d * sinb * sin(phi_d) - z = r_d * cosb - v_d = np.dstack((x, y, z)) - - obj = cls(v_d, id="pxrdip") - obj.xs = xs - obj.twotheta_n = twotheta_n - obj.phi = phi - obj.centered_ph = centered_ph - return obj - - @classmethod - def from_point_params(cls, r_d=24, beta=45, phi_d=0, id=None): - """Single-point detector given by spherical coordinates.""" - beta = radians(beta) - phi_d = radians(phi_d) - v_d = r_d * np.array((sin(beta)*cos(phi_d), sin(beta)*sin(phi_d), cos(beta))) - return cls(v_d, id) - - @property - def shape(self): - return self.beta.shape - - def calc_twotheta_n(self, xs): - """Return nominal twotheta for all elements for given XS and phi=0 vector.""" - return calc_twotheta_nominal(xs.alpha, self.beta, xs.phi, self.phi_d) - - def calc_phi(self, xs, v_phi0=(0,0,1)): - """Return phi for all elements for given XS and phi=0 vector.""" - v_phi0 = np.array(v_phi0) - v_x = xs.v_x - v_d = self.v_d - - phi = azimuth(v_d, v0=v_x, v1=v_phi0) - return phi - - -class ReportFileWriter(): - """Writes out the twotheta_pinhole correction report file with tables.""" - TWOTHETA_MIN = 10 # [deg] no useful data this close to direct image - BETA_MIN = 10 # [deg] VISAR hole is 25 mm diam at 75 mm distance =~10 deg - BETA_MAX = 88 # [deg] correction value very close to 90 deg is dubious - DEFAULT_NP = 240 # default number of pinhole phi points for integration - - def __init__(self, params=None): - p = params if isinstance(params, dict) else EP_BASELINE - self.params = p - self.xs = XraySource(p['mat_x'], p['r_x'], p['alpha'], p['phi_x']) - self.ph = Pinhole(p['mat_p'], p['d_p'], p['h_p']) - self.centered_ph = p['centered_ph'] - - self.beta_min, self.beta_max = self.BETA_MIN, self.BETA_MAX - self.Np = self.DEFAULT_NP - self.mu_p = MU_P_LOOKUP.get((self.ph.mat, self.xs.mat), 1/3.0) # [um^-1] attenuation coefficient - self.run_calc() - - @property - def title(self): - xs, ph = self.xs, self.ph - return "{}{}_a{:.3g}_r{:.0f}_d{:.0f}_h{:.0f}".format(xs.mat, ph.mat, xs.alpha_deg, xs.r, 1000*ph.diam, 1000*ph.thick) - - def run_calc(self): - if True: # standard 2theta and phi sample spacing - d2theta, dphi = 2, 5 - else: # extra fine 2theta and phi sample spacing - d2theta, dphi = 0.5, 1 - - # create detector, and generate range for twotheta and phi - xs, ph = self.xs, self.ph - qq = np.arange(self.TWOTHETA_MIN, xs.alpha_deg + self.beta_max, d2theta) - pp = np.arange(0, 180 + 0.001, dphi) - twotheta_n, phi = np.meshgrid(qq, pp) - det = Detector.pxrdip_from_twothetaphi(xs, twotheta_n, phi, self.centered_ph) - - self.qq, self.pp = qq, pp - self.twotheta_n, self.phi = twotheta_n, phi - self.det = det - - qq_n = det.twotheta_n - qq_p, V_p = calc_twotheta_pinhole(xs, ph, det, Np=self.Np) - qq_JHE = calc_twotheta_JHE(xs, ph, det) - - qq_p[det.beta < np.radians(self.beta_min)] = np.nan - qq_p[det.beta > np.radians(self.beta_max)] = np.nan - dq_pn = 180 / pi * (qq_p - qq_n) # convert to degrees - dq_pj = 180 / pi * (qq_p - qq_JHE) # convert to degrees - - self.qq_n, self.qq_p = qq_n, qq_p - self.dq_pn, self.dq_pj = dq_pn, dq_pj - - if False: # switch on/off additional tables - write_table(dq_pn.T, self.title + "_pn") - write_table(dq_pj.T, self.title + "_pj") - write_table(180/pi*qq_n[0,:][None,:], self.title + "_qq_n", fmt="%.0f") - write_table(180/pi*det.phi[:,0][None,:], self.title + "_phi", fmt="%.0f") - - - def write(self, filename=None, tables="standard"): - """Write pxrdip twotheta pinhole correction tables to file""" - if filename is None: - name = "pxrdip-2theta-pinhole-correction_{}.txt".format(self.title) - filename = osp.join(OUTPUT_PATH, name) - if tables == "standard": - self.tables = tables = ( - (self.dq_pn.T, "%6.3f", "twotheta_p - twotheta_n", "deg"), - (self.dq_pj.T, "%6.3f", "twotheta_p - twotheta_JHE", "deg"), - ) - ss = [HEADER] - ss.append("[parameters]") - ss.extend(self.param_output_list()) - - ss.append("\n[{} rows, mapped to twotheta_n, deg]".format(len(self.qq))) - ss.append(" ".join(["{:.1f}".format(q) for q in self.qq])) - ss.append("\n[{} columns, mapped to phi, deg]".format(len(self.pp))) - ss.append(" ".join(["{:.0f}".format(q) for q in self.pp])) - - for i, table in enumerate(tables): - array, fmt, label, uom = table - ss.append("\n[table{}, {}, {}]".format(i, label, uom)) - with io.StringIO() as sio: - np.savetxt(sio, array, fmt=fmt) - ss.append(sio.getvalue()) - - with open(filename, "w") as f: - f.write("\n".join(ss)) - - def param_output_list(self): - d = self.params - ss = [] - ss.append("xsf: mat={mat_x}, r_x={r_x} mm, alpha={alpha} deg, phi_x={phi_x} deg".format(**d)) - ss.append("pinhole: mat={}, d_p={:.0f} um, h_p={:.0f} um".format(d['mat_p'], 1000*d['d_p'], 1000*d['h_p'])) - ss.append("attenuation coefficient: mu_p=1/({:.2f} um)".format(1/self.mu_p)) - ss.append("diagnostic: pxrdip, centered pinhole") - ss.append("integration_resolution: Np={}".format(self.Np)) - ss.append("evaluation_script: {} version {}".format(osp.split(__file__)[1], __version__)) - ss.append("evaluation_timestamp: {}".format(datetime.datetime.today())) - ss.append("row, column axes: twotheta_n, phi") - - for i, table in enumerate(self.tables): - array, fmt, label, uom = table - ss.append("table{}: {}".format(i, label)) - return ss - - def save_image(self, filename=None): - if filename is None: - name = "pxrdip-2theta-pinhole-correction_{}.png".format(self.title) - filename = osp.join(OUTPUT_PATH, name) - - xs, ph, det = self.xs, self.ph, self.det - qq_n = self.qq_n - dq_pn, dq_pj = self.dq_pn, self.dq_pj - - extent = 180/pi*np.array((np.nanmin(qq_n), np.nanmax(qq_n), np.min(det.phi), np.max(det.phi))) - extent2 = extent * np.array((1,1,-1,-1)) - kws = dict(extent=extent, aspect='auto', cmap='coolwarm') - kws2 = dict(extent=extent2, aspect='auto', cmap='coolwarm') - kws_contour = dict(colors='k', linestyles='--') - - def show_image(ax, value, label): - vmax = max(np.nanmax(value), -np.nanmin(value)) - im = ax.imshow(value, vmin=-vmax, vmax=vmax, **kws2) - im = ax.imshow(value, vmin=-vmax, vmax=vmax, **kws) - plt.colorbar(im, ax=ax, label=label) - ax.contour(det.beta, [ph.critical_angle], extent=extent, **kws_contour) - ax.contour(det.beta, [ph.critical_angle], extent=extent2, **kws_contour) - - thetalabel, philabel = r'$2\theta_n$ [deg]', r'$\phi$ [deg]' - fig = plt.figure(figsize=(5,6), dpi=DPI_SCREEN) - ax1 = plt.subplot(211, xlabel=thetalabel, ylabel=philabel, title=" PXRDIP {}".format(self.title)) - ax2 = plt.subplot(212, sharex=ax1, sharey=ax1, xlabel=thetalabel, ylabel=philabel) - show_image(ax1, dq_pn, r'$2\theta_p$ - $2\theta_n$ [deg]') - show_image(ax2, dq_pj, r'$2\theta_p$ - $2\theta_{JHE}$ [deg]') - - # formatting tweaks - ax1.set_xticks(np.arange(0,151,15)) - ax1.set_xticks(np.arange(0,151,5), minor=True) - ax1.set_yticks(np.arange(-180,181,90)) - ax1.set_yticks(np.arange(-180,181,45), minor=True) - ax1.axis((self.TWOTHETA_MIN, xs.alpha_deg + self.beta_max, -180, 180)) - plt.tight_layout() - - if False: # optional: toggle printing additional info - print("shape (Nphi, N2theta) and size of output array:", np.shape(dq_pj), np.size(dq_pj)) - print("max difference in twotheta direction: {:.4f} deg".format(np.nanmax(dq_pj[:,1:]-dq_pj[:,:-1]))) - print("max difference in phi direction: {:.4f} deg".format(np.nanmax(dq_pj[1:,:]-dq_pj[:-1,:]))) - print("max correction {:.4f} deg".format(np.nanmax(dq_pj[1:,:]))) - print("min correction {:.4f} deg".format(np.nanmin(dq_pj[1:,:]))) - - fig.savefig(filename, dpi=200) - - -def write_table(A, filename=None, fmt="%6.3f"): - """write table to file""" - if filename is None: - filename = osp.join(OUTPUT_PATH, "output.txt") - else: - filename = osp.join(OUTPUT_PATH, filename + ".txt") - np.savetxt(filename, A, fmt=fmt) - - -def sandbox1(): - """Sandbox for exploring twotheta calculations.""" - # speed vs fidelity parameters - calc_type = 1 - if calc_type in ("fast", 0): - d2theta, dphi = 5, 15 # 2theta and phi sample spacing - Np = 72 # number of azimuthal points around pinhole - elif calc_type in ("hires", 2): - d2theta, dphi = 0.5, 1 # 2theta and phi sample spacing - Np = 240 # number of azimuthal points around pinhole - else: - d2theta, dphi = 2, 5 # 2theta and phi sample spacing - Np = 240 # number of azimuthal points around pinhole - - twothetamin = 10 # [deg] no useful data this close to direct image - betamin, betamax = 10, 88 # deg - - # xray source and pinhole - if False: -# p = OMEGA_BASELINE - p = EP_BASELINE - else: -# p = dict(mat_x='Cu', r_x=24.24, alpha=22.5, phi_x=0, -# mat_p='Ta', d_p=0.3, h_p=0.075, centered_ph=True) -# p = dict(mat_x='Cu', r_x=24.24, alpha=22.5, phi_x=0, z_s=0.1, -# mat_p='Ta', d_p=0.8, h_p=0.075, centered_ph=True) # Danae sodium, although actually with non-centered pinhole - p = dict(mat_x='Cu', r_x=17, alpha=45, phi_x=0, z_s=0.1, - mat_p='W', d_p=0.5, h_p=0.075, centered_ph=True) # Michelle TATB - - xs = XraySource(p['mat_x'], p['r_x'], p['alpha'], p['phi_x']) - ph = Pinhole(p['mat_p'], p['d_p'], p['h_p']) - - # create detector, and generate range for twotheta and phi - qq = np.arange(twothetamin, xs.alpha_deg + betamax, d2theta) - pp = np.arange(0, 180 + 1, dphi) - twotheta_n, phi = np.meshgrid(qq, pp) - det = Detector.pxrdip_from_twothetaphi(xs, twotheta_n, phi, p['centered_ph']) - - title = "{}{}_a{:.3g}_d{:.0f}".format(xs.mat, ph.mat, xs.alpha_deg, 1000*ph.diam) - - qq_n = det.twotheta_n - qq_p, V_p = calc_twotheta_pinhole(xs, ph, det, Np=Np) - qq_JHE = calc_twotheta_JHE(xs, ph, det) - - qq_p[det.betanp.radians(betamax)] = np.nan - dq_pn = 180/pi*(qq_p - qq_n) - dq_pj = 180/pi*(qq_p - qq_JHE) - - extent = 180/pi*np.array((np.nanmin(qq_n), np.nanmax(qq_n), np.min(det.phi), np.max(det.phi))) - extent2 = extent * np.array((1,1,-1,-1)) - kws = dict(extent=extent, aspect='auto', cmap='coolwarm') - kws2 = dict(extent=extent2, aspect='auto', cmap='coolwarm') - kws_contour = dict(colors='k', linestyles='--') - - - def show_image(ax, value, label): - vmax = max(np.nanmax(value), -np.nanmin(value)) - im = ax.imshow(value, vmin=-vmax, vmax=vmax, **kws2) - im = ax.imshow(value, vmin=-vmax, vmax=vmax, **kws) - plt.colorbar(im, ax=ax, label=label) - ax.contour(det.beta, [ph.critical_angle], extent=extent, **kws_contour) - ax.contour(det.beta, [ph.critical_angle], extent=extent2, **kws_contour) - - plt.figure(figsize=(6,7), dpi=DPI_SCREEN) - ax1 = plt.subplot(211, xlabel=r'$2\theta_n$', ylabel=r'$\phi$', title=title) - ax2 = plt.subplot(212, sharex=ax1, sharey=ax1, xlabel=r'$2\theta_n$', ylabel=r'$\phi$') - show_image(ax1, dq_pn, r'$2\theta_p$ - $2\theta_n$') - show_image(ax2, dq_pj, r'$2\theta_p$ - $2\theta_{JHE}$') - - # formatting tweaks -# [ax.minorticks_on() for ax in (ax1, ax2)] - ax1.set_xticks(np.arange(0,151,15)) - ax1.set_xticks(np.arange(0,151,5), minor=True) - ax1.set_yticks(np.arange(-180,181,90)) - ax1.set_yticks(np.arange(-180,181,45), minor=True) - ax1.axis((twothetamin, xs.alpha_deg + betamax, -180, 180)) - plt.tight_layout() - - print("shape (Nphi, N2theta) and size of output array:", np.shape(dq_pj), np.size(dq_pj)) - print("max difference in twotheta direction: {:.4f} deg".format(np.nanmax(dq_pj[:,1:]-dq_pj[:,:-1]))) - print("max difference in phi direction: {:.4f} deg".format(np.nanmax(dq_pj[1:,:]-dq_pj[:-1,:]))) - - if calc_type not in ("hires", 2): - write_table(dq_pn.T, title + "_pn") - write_table(dq_pj.T, title + "_pj") - write_table(180/pi*qq_n[0,:][None,:], title + "_qq_n", fmt="%.0f") - write_table(180/pi*det.phi[:,0][None,:], title + "_phi", fmt="%.0f") - - if True: # impact parameter calc - z_s = p['z_s'] - b_center = z_s * np.tan(det.beta) - b_max = ph.d/2 + (z_s - ph.h/2) * np.tan(det.beta) - b_center[det.beta >= ph.critical_angle] = np.nan - b_max[det.beta >= ph.critical_angle] = np.nan - plt.figure(figsize=(6,7), dpi=DPI_SCREEN) - ax1 = plt.subplot(211, xlabel=r'$2\theta_n$', ylabel=r'$\phi$', title=title) - ax2 = plt.subplot(212, sharex=ax1, sharey=ax1, xlabel=r'$2\theta_n$', ylabel=r'$\phi$') - - im = ax1.imshow(b_center, vmin=0, **kws2) - im = ax1.imshow(b_center, vmin=0, **kws) - plt.colorbar(im, ax=ax1, label="b center [mm]") - im = ax2.imshow(b_max, vmin=ph.d/2, **kws2) - im = ax2.imshow(b_max, vmin=ph.d/2, **kws) - plt.colorbar(im, ax=ax2, label="b max [mm]") - - ax1.contour(det.beta, [ph.critical_angle], extent=extent, **kws_contour) - ax1.contour(det.beta, [ph.critical_angle], extent=extent2, **kws_contour) - ax1.set_xticks(np.arange(0,151,15)) - ax1.set_xticks(np.arange(0,151,5), minor=True) - ax1.set_yticks(np.arange(-180,181,90)) - ax1.set_yticks(np.arange(-180,181,45), minor=True) - ax1.axis((twothetamin, xs.alpha_deg + betamax, -180, 180)) - - plt.tight_layout() - - if calc_type not in ("hires", 2): - write_table(b_center.T, title + "_b-center") - if calc_type not in ("hires", 2): - write_table(b_max.T, title + "_b-max") - -# ============================ MAIN =========================================== -if __name__ == "__main__": - if True: - p = OMEGA_BASELINE - else: - p = EP_BASELINE - rfw = ReportFileWriter(p) - rfw.write() - rfw.save_image() - -# sandbox1() - - plt.show() diff --git a/scripts/plot_instrument_schematic.py b/scripts/plot_instrument_schematic.py deleted file mode 100644 index 9e5c51558..000000000 --- a/scripts/plot_instrument_schematic.py +++ /dev/null @@ -1,179 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Wed Oct 26 12:10:16 2022 - -@author: jbernier -""" - -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Thu Aug 4 16:18:20 2022 - -@author: jbernier -""" -import os - -import numpy as np - -from mpl_toolkits.mplot3d import Axes3D -from mpl_toolkits.mplot3d.proj3d import proj_transform -from mpl_toolkits.mplot3d.art3d import Poly3DCollection - -from matplotlib.text import Annotation -from matplotlib.patches import FancyArrowPatch - -import matplotlib.pyplot as plt - -import yaml - -from hexrd import instrument - -# %% -hexrd_dir = '/home/jbernier/Documents/GitHub/hexrd' -icfg = yaml.safe_load( - open(os.path.join(hexrd_dir, - 'hexrd/resources/pxrdip_reference_config.yml'), - 'r') - ) - -instr = instrument.HEDMInstrument(icfg) -instr.source_distance = 22.5 - -colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] - -# %% -class Annotation3D(Annotation): - - def __init__(self, text, xyz, *args, **kwargs): - super().__init__(text, xy=(0, 0), *args, **kwargs) - self._xyz = xyz - - def draw(self, renderer): - x2, y2, z2 = proj_transform(*self._xyz, self.axes.M) - self.xy = (x2, y2) - super().draw(renderer) - -class Arrow3D(FancyArrowPatch): - - def __init__(self, x, y, z, dx, dy, dz, *args, **kwargs): - super().__init__((0, 0), (0, 0), *args, **kwargs) - self._xyz = (x, y, z) - self._dxdydz = (dx, dy, dz) - - def draw(self, renderer): - x1, y1, z1 = self._xyz - dx, dy, dz = self._dxdydz - x2, y2, z2 = (x1 + dx, y1 + dy, z1 + dz) - - xs, ys, zs = proj_transform((x1, x2), (y1, y2), (z1, z2), self.axes.M) - self.set_positions((xs[0], ys[0]), (xs[1], ys[1])) - super().draw(renderer) - - def do_3d_projection(self, renderer=None): - x1, y1, z1 = self._xyz - dx, dy, dz = self._dxdydz - x2, y2, z2 = (x1 + dx, y1 + dy, z1 + dz) - - xs, ys, zs = proj_transform((x1, x2), (y1, y2), (z1, z2), self.axes.M) - self.set_positions((xs[0], ys[0]), (xs[1], ys[1])) - - return np.min(zs) - - -def _annotate3D(ax, text, xyz, *args, **kwargs): - '''Add anotation `text` to an `Axes3d` instance.''' - - annotation = Annotation3D(text, xyz, *args, **kwargs) - ax.add_artist(annotation) - - -setattr(Axes3D, 'annotate3D', _annotate3D) - - -def _arrow3D(ax, x, y, z, dx, dy, dz, *args, **kwargs): - '''Add an 3d arrow to an `Axes3D` instance.''' - - arrow = Arrow3D(x, y, z, dx, dy, dz, *args, **kwargs) - ax.add_artist(arrow) - - -setattr(Axes3D, 'arrow3D', _arrow3D) - -# %% -fig = plt.figure() -# ax = Axes3D(fig, auto_add_to_figure=False) -# fig.add_axes(ax) -ax = fig.add_subplot(111, projection='3d') - -rho = instr.source_distance - - -def draw_cs(ax, origin=np.zeros(3), scale=10, mscale=10, rmat=None): - dirs = np.hstack([np.tile(origin.flatten(), (3, 1)), - [i for i in scale*np.eye(3)]]) - if rmat is not None: - for i in range(3): - dirs[i, -3:] = np.dot(dirs[i, -3:], rmat.T) - - colors = ['r', 'g', 'b'] - for i, qargs in enumerate(dirs): - # ax.quiver3D(*qargs, color=colors[i], linewidth=2) - ax.arrow3D(*qargs, - mutation_scale=mscale, - ec=colors[i], - fc=colors[i]) - - - -for i, (det_name, panel) in enumerate(instr.detectors.items()): - verts = np.dot( - np.hstack( - [np.vstack([panel.corner_ll, - panel.corner_lr, - panel.corner_ur, - panel.corner_ul, - panel.corner_ll]), - np.zeros((5, 1))] - ), - panel.rmat.T) + panel.tvec - ax.add_collection3d( - Poly3DCollection( - [[tuple(i) for i in verts]], - facecolor=colors[i], edgecolor='k', alpha=0.5 - ), zdir='y' - ) - #ax.quiver3D(0, 0, 0, *panel.tvec, color='m') - ax.arrow3D(0, 0, 0, *panel.tvec, - mutation_scale=10, - ec ='m', - fc='c') - ax.annotate3D(det_name, panel.tvec, - xytext=(3, 3), textcoords='offset points', color=colors[i]) - draw_cs(ax, panel.tvec, rmat=panel.rmat, mscale=5, scale=15) - -draw_cs(ax, origin=np.zeros(3), mscale=20, scale=50) - - -# XRS plotting -xrs_vec = -rho*instr.beam_vector -# ax.quiver3D(0, 0, 0, *xrs_vec, color='k') -ax.arrow3D(0, 0, 0, *xrs_vec, - mutation_scale=10, - ec ='m', - fc='k') -ax.scatter3D(*xrs_vec, marker='o', color='k', s=48) - -# scene -ax.set_xlim(-100, 100) -ax.set_ylim(-100, 100) -ax.set_zlim(-100, 100) -ax.set_xlabel(r'$\mathbf{X}_l$', fontsize=24, color='r') -ax.set_ylabel(r'$\mathbf{Y}_l$', fontsize=24, color='g') -ax.set_zlabel(r'$\mathbf{Z}_l$', fontsize=24, color='b') - -fig.tight_layout() - -plt.draw() -plt.show() diff --git a/scripts/powder_calibration.py b/scripts/powder_calibration.py deleted file mode 100644 index 5fac2db3f..000000000 --- a/scripts/powder_calibration.py +++ /dev/null @@ -1,249 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Jul 30 20:45:01 2019 - -@author: joel -""" - -import numpy as np - -from scipy.optimize import least_squares, leastsq - -from hexrd.fitting import fitpeak -from hexrd.matrixutil import findDuplicateVectors - - -# %% -class InstrumentCalibrator(object): - def __init__(self, *args): - assert len(args) > 0, \ - "must have at least one calibrator" - self._calibrators = args - self._instr = self._calibrators[0].instr - - @property - def instr(self): - return self._instr - - @property - def calibrators(self): - return self._calibrators - - # ========================================================================= - # METHODS - # ========================================================================= - - def run_calibration(self, use_robust_optimization=False): - """ - FIXME: only coding serial powder case to get things going. Will - eventually figure out how to loop over multiple calibrator classes. - All will have a reference the same instrument, but some -- like single - crystal -- will have to add parameters as well as contribute to the RHS - """ - calib_class = self.calibrators[0] - - obj_func = calib_class.residual - - data_dict = calib_class._extract_powder_lines() - - # grab reduced optimizaion parameter set - x0 = self._instr.calibration_parameters[ - self._instr.calibration_flags - ] - - resd0 = obj_func(x0, data_dict) - - if use_robust_optimization: - oresult = least_squares( - obj_func, x0, args=(data_dict, ), - method='trf', loss='soft_l1' - ) - x1 = oresult['x'] - else: - x1, cox_x, infodict, mesg, ierr = leastsq( - obj_func, x0, args=(data_dict, ), - full_output=True - ) - resd1 = obj_func(x1, data_dict) - - delta_r = sum(resd0**2)/float(len(resd0)) - \ - sum(resd1**2)/float(len(resd1)) - - if delta_r > 0: - print(('OPTIMIZATION SUCCESSFUL\nfinal ssr: %f' % sum(resd1**2))) - print(('delta_r: %f' % delta_r)) - # self.instr.write_config(instrument_filename) - else: - print('no improvement in residual!!!') - - -# %% -class PowderCalibrator(object): - def __init__(self, instr, plane_data, img_dict, - tth_tol=None, eta_tol=0.25, - pktype='pvoigt'): - assert list(instr.detectors.keys()) == list(img_dict.keys()), \ - "instrument and image dict must have the same keys" - self._instr = instr - self._plane_data = plane_data - self._img_dict = img_dict - - # for polar interpolation - self._tth_tol = tth_tol or np.degrees(plane_data.tThWidth) - self._eta_tol = eta_tol - - # for peak fitting - # ??? fitting only, or do alternative peak detection? - self._pktype = pktype - - @property - def instr(self): - return self._instr - - @property - def plane_data(self): - return self._plane_data - - @property - def img_dict(self): - return self._img_dict - - @property - def tth_tol(self): - return self._tth_tol - - @tth_tol.setter - def tth_tol(self, x): - assert np.isscalar(x), "tth_tol must be a scalar value" - self._tth_tol = x - - @property - def eta_tol(self): - return self._eta_tol - - @eta_tol.setter - def eta_tol(self, x): - assert np.isscalar(x), "eta_tol must be a scalar value" - self._eta_tol = x - - @property - def pktype(self): - return self._pktype - - @pktype.setter - def pktype(self, x): - """ - currently only 'pvoigt' or 'gaussian' - """ - assert isinstance(x, str), "tth_tol must be a scalar value" - self._pktype = x - - def _interpolate_images(self): - """ - returns the iterpolated powder line data from the images in img_dict - - ??? interpolation necessary? - """ - return self.instr.extract_line_positions( - self.plane_data, self.img_dict, - tth_tol=self.tth_tol, eta_tol=self.eta_tol, - npdiv=2, collapse_eta=False, collapse_tth=False, - do_interpolation=True) - - def _extract_powder_lines(self): - """ - return the RHS for the instrument DOF and image dict - - The format is a dict over detectors, each containing - - [index over ring sets] - [index over azimuthal patch] - [xy_meas, tth_meas, tth_ref, eta_ref] - - FIXME: can not yet handle tth ranges with multiple peaks! - """ - # ideal tth - tth_ideal = self.plane_data.getTTh() - tth0 = [] - for idx in self.plane_data.getMergedRanges()[0]: - if len(idx) > 1: - eqv, uidx = findDuplicateVectors(np.atleast_2d(tth_ideal[idx])) - if len(uidx) > 1: - raise NotImplementedError("can not handle multipeak yet") - else: - # if here, only degenerate ring case - uidx = idx[0] - else: - uidx = idx[0] - tth0.append(tth_ideal[uidx]) - - powder_lines = self._interpolate_images() - - # GRAND LOOP OVER PATCHES - rhs = dict.fromkeys(self.instr.detectors) - for det_key, panel in self.instr.detectors.items(): - rhs[det_key] = [] - for i_ring, ringset in enumerate(powder_lines[det_key]): - tmp = [] - for angs, intensities in ringset: - tth_centers = np.average( - np.vstack([angs[0][:-1], angs[0][1:]]), - axis=0) - eta_ref = angs[1] - int1d = np.sum(np.array(intensities).squeeze(), axis=0) - """ - DARREN: FIT [tth_centers, intensities[0]] HERE - - RETURN TTH0 - rhs.append([tth0, eta_ref]) - """ - p0 = fitpeak.estimate_pk_parms_1d( - tth_centers, int1d, self.pktype - ) - - p = fitpeak.fit_pk_parms_1d( - p0, tth_centers, int1d, self.pktype - ) - # !!! this is where we can kick out bunk fits - tth_meas = p[1] - center_err = abs(tth_meas - tth0[i_ring]) - if p[0] < 0.1 or center_err > np.radians(self.tth_tol): - continue - xy_meas = panel.angles_to_cart([[tth_meas, eta_ref], ]) - tmp.append( - np.hstack( - [xy_meas.squeeze(), - tth_meas, - tth0[i_ring], - eta_ref] - ) - ) - rhs[det_key].append(np.vstack(tmp)) - rhs[det_key] = np.vstack(rhs[det_key]) - return rhs - - def residual(self, reduced_params, data_dict): - """ - """ - - # first update instrument from input parameters - full_params = self.instr.calibration_parameters - full_params[self.instr.calibration_flags] = reduced_params - self.instr.update_from_parameter_list(full_params) - - # build residual - resd = [] - for det_key, panel in self.instr.detectors.items(): - pdata = np.vstack(data_dict[det_key]) - if len(pdata) > 0: - calc_xy = panel.angles_to_cart(pdata[:, -2:]) - - resd.append( - (pdata[:, :2].flatten() - calc_xy.flatten()) - ) - else: - continue - return np.hstack(resd) - - diff --git a/scripts/preprocess_tomo_mask.py b/scripts/preprocess_tomo_mask.py deleted file mode 100644 index 8f2bcdb81..000000000 --- a/scripts/preprocess_tomo_mask.py +++ /dev/null @@ -1,214 +0,0 @@ -#%% Necessary Dependencies - -import numpy as np -import logging - -import yaml - -try: - import matplotlib.pyplot as plt - matplot=True -except(ImportError): - logging.warning(f'no matplotlib, debug plotting disabled') - matplot=False - -from hexrd.grainmap import nfutil -from hexrd.grainmap import tomoutil - - -from hexrd import instrument - - -def load_instrument(yml): - with open(yml, 'r') as f: - icfg = yaml.load(f, Loader=yaml.FullLoader) - return instrument.HEDMInstrument(instrument_config=icfg) - -# %% FILES TO LOAD -CAN BE EDITED -#============================================================================== -#These files are attached, retiga.yml is a detector configuration file -#The near field detector was already calibrated - -#A materials file, is a cPickle file which contains material information like lattice -#parameters necessary for the reconstruction - -main_dir = '/INSERT/WORK/DIR/' - -det_file = main_dir + 'tomo_det.yml' - -#============================================================================== -# %% OUTPUT INFO -CAN BE EDITED -#============================================================================== - -output_dir = main_dir -output_stem='tomo_out' - -#============================================================================== -# %% TOMOGRAPHY DATA FILES -CAN BE EDITED -#============================================================================== - -stem='nf_' - -#Locations of tomography dark field images -tdf_data_folder='/LOC/nf/' - -tdf_img_start=52 #for this rate, this is the 6th file in the folder -tdf_num_imgs=10 - -#Locations of tomography bright field images -tbf_data_folder='/LOC/nf/' - -tbf_img_start=68 #for this rate, this is the 6th file in the folder -tbf_num_imgs=10 - -#Locations of tomography images -tomo_data_folder='/LOC/nf/' - -tomo_img_start=84#for this rate, this is the 6th file in the folder -tomo_num_imgs=360 - -#============================================================================== -# %% USER OPTIONS -CAN BE EDITED -#============================================================================== - - -ome_range_deg=[(0.,359.75)] #degrees - - -#tomography options -recon_thresh=0.0002#usually varies between 0.0001 and 0.0005 -#Don't change these unless you know what you are doing, this will close small holes -#and remove noise -noise_obj_size=500 -min_hole_size=500 -erosion_iter=1 -dilation_iter=1 -project_single_layer=False #projects the center layers through the volume, faster but not recommended, included for completion / historical purposes - - -#reconstruction volume options -cross_sectional_dim=1.35 #cross sectional to reconstruct (should be at least 20%-30% over sample width) -voxel_spacing=0.005#in mm -v_bnds=[-0.4,0.4] - - -#============================================================================== -# %% LOAD INSTRUMENT DATA -#============================================================================== - - -instr=load_instrument(det_file) - -panel = next(iter(instr.detectors.values())) - - - -nrows=panel.rows -ncols=panel.cols -pixel_size=panel.pixel_size_row - - -rot_axis_pos=panel.tvec[0] #should match t_vec_d[0] from nf_detector_parameter_file -vert_beam_center=panel.tvec[1] - - -# need to do a few calculations because not every row will be reconstructed -# depending on sampling -vert_points=np.arange(v_bnds[0]+voxel_spacing/2.,v_bnds[1],voxel_spacing) - -center_layer_row=nrows/2.+vert_beam_center/pixel_size - -rows_to_recon=np.round(center_layer_row-vert_points/pixel_size).astype(int) - -center_layer_row=int(center_layer_row) - -#============================================================================== -# %% TOMO PROCESSING - GENERATE DARK AND BRIGHT FIELD -#============================================================================== -tdf=tomoutil.gen_median_image(tdf_data_folder,tdf_img_start,tdf_num_imgs,nrows,ncols,stem=stem,num_digits=6) -tbf=tomoutil.gen_median_image(tbf_data_folder,tbf_img_start,tbf_num_imgs,nrows,ncols,stem=stem,num_digits=6) - -#============================================================================== -# %% TOMO PROCESSING - BUILD RADIOGRAPHS -#============================================================================== - -rad_stack=tomoutil.gen_attenuation_rads(tomo_data_folder,tbf,tomo_img_start,tomo_num_imgs,nrows,ncols,stem=stem,num_digits=6,tdf=tdf) - - - -#============================================================================== -# %% TOMO PROCESSING - INVERT SINOGRAM -#============================================================================== -# center = 0.0 - -test_fbp=tomoutil.tomo_reconstruct_layer(rad_stack,cross_sectional_dim,layer_row=center_layer_row,\ - start_tomo_ang=ome_range_deg[0][0],end_tomo_ang=ome_range_deg[0][1],\ - tomo_num_imgs=tomo_num_imgs, center=rot_axis_pos,pixel_size=pixel_size) - -test_binary_recon=tomoutil.threshold_and_clean_tomo_layer(test_fbp,recon_thresh, \ - noise_obj_size,min_hole_size, erosion_iter=erosion_iter, \ - dilation_iter=dilation_iter) - -tomo_mask_center=tomoutil.crop_and_rebin_tomo_layer(test_binary_recon,recon_thresh,voxel_spacing,pixel_size,cross_sectional_dim) - - -#============================================================================== -# %% TOMO PROCESSING - VIEW RAW FILTERED BACK PROJECTION -#============================================================================== - -if matplot: - plt.figure(1) - plt.imshow(test_fbp,vmin=recon_thresh,vmax=recon_thresh*2) - plt.title('Check Thresholding') - #Use this image to view the raw reconstruction, estimate threshold levels. and - #figure out if the rotation axis position needs to be corrected - - - plt.figure(2) - plt.imshow(tomo_mask_center,interpolation='none') - plt.title('Check Center Mask') - -#============================================================================== -# %% PROCESS REMAINING LAYERS -#============================================================================== - -full_mask=np.zeros([len(rows_to_recon),tomo_mask_center.shape[0],tomo_mask_center.shape[1]]) - - -for ii in np.arange(len(rows_to_recon)): - print('Layer: ' + str(ii) + ' of ' + str(len(rows_to_recon))) - - if project_single_layer: #not recommended option - full_mask[ii]=tomo_mask_center - - else: - reconstruction_fbp=tomoutil.tomo_reconstruct_layer(rad_stack,cross_sectional_dim,layer_row=rows_to_recon[ii],\ - start_tomo_ang=ome_range_deg[0][0],end_tomo_ang=ome_range_deg[0][1],\ - tomo_num_imgs=tomo_num_imgs, center=rot_axis_pos,pixel_size=pixel_size) - - binary_recon=tomoutil.threshold_and_clean_tomo_layer(reconstruction_fbp,recon_thresh, \ - noise_obj_size,min_hole_size,erosion_iter=erosion_iter, \ - dilation_iter=dilation_iter) - - tomo_mask=tomoutil.crop_and_rebin_tomo_layer(binary_recon,recon_thresh,voxel_spacing,pixel_size,cross_sectional_dim) - - full_mask[ii]=tomo_mask - -#============================================================================== -# %% TOMO PROCESSING - VIEW LAST TOMO_MASK FOR SAMPLE BOUNDS -#============================================================================== - -if matplot: - plt.figure(3) - plt.imshow(tomo_mask,interpolation='none') - plt.title('Check Center Mask') - -#============================================================================== -# %% TOMO PROCESSING - CONSTRUCT DATA GRID -#============================================================================== - -test_crds, n_crds, Xs, Ys, Zs = nfutil.gen_nf_test_grid_tomo(full_mask.shape[2], full_mask.shape[1], v_bnds, voxel_spacing) - -#%% - -np.savez('tomo_mask.npz',mask=full_mask,Xs=Xs,Ys=Ys,Zs=Zs,voxel_spacing=voxel_spacing) diff --git a/scripts/process_nf.py b/scripts/process_nf.py deleted file mode 100644 index 41dcf35bb..000000000 --- a/scripts/process_nf.py +++ /dev/null @@ -1,364 +0,0 @@ -#!/usr/bin/env python2 -# -*- coding: utf-8 -*- -""" - -@author: dcp5303 -""" - -#%% Necessary Dependencies - - -# PROCESSING NF GRAINS WITH MISORIENTATION -#============================================================================== -import numpy as np - -import multiprocessing as mp - -import os -import logging -import psutil - -from hexrd import rotations as rot -from hexrd.grainmap import nfutil -from hexrd.material import Material -from hexrd.valunits import valWUnit - -try: - import matplotlib.pyplot as plt - matplot = True -except(ImportError): - logging.warning(f'no matplotlib, debug plotting disabled') - matplot = False - - -#============================================================================== -# %% FILES TO LOAD -CAN BE EDITED -#============================================================================== -#These files are attached, retiga.yml is a detector configuration file -#The near field detector was already calibrated - -#A materials file, is a cPickle file which contains material information like lattice -#parameters necessary for the reconstruction - - -main_dir = '/INSERT/WORKDIR/' - -det_file = main_dir + 'retiga.yml' -mat_file = main_dir + 'materials.h5' - -#============================================================================== -# %% OUTPUT INFO -CAN BE EDITED -#============================================================================== - -output_dir = main_dir - - -#============================================================================== -# %% NEAR FIELD DATA FILES -CAN BE EDITED - ZERO LOAD SCAN -#============================================================================== -#These are the near field data files used for the reconstruction, a grains.out file -#from the far field analaysis is used as orientation guess for the grid that will -grain_out_file = '/LOC/grains.out' - -#%% - -stem = 'nf_' - -#Locations of near field images -data_folder = '/INSERT/nf/' - -img_start = 6 # whatever you want the first frame to be -nframes = 1440 -img_nums = np.arange(img_start, img_start+nframes, 1) - -output_stem = 'NAME_VOL_X' -#============================================================================== -# %% USER OPTIONS -CAN BE EDITED -#============================================================================== - -beam_energy = 71.667 - - -### material for the reconstruction -mat_name = 'INSERT' -max_tth = None # degrees, if None is input max tth will be set by the geometry - - -#reconstruction with misorientation included, for many grains, this will quickly -#make the reconstruction size unmanagable -misorientation_bnd = 0.0 # degrees -misorientation_spacing = 0.25 # degrees - - -#####image processing parameters -num_for_dark = 250 # num images to use for median data - -img_threshold = 0 -process_type = 'gaussian' -sigma = 2.0 -size = 3.0 -process_args = np.array([sigma, size]) - -# process_type='dilations_only' -# num_erosions=2 #num iterations of images erosion, don't mess with unless you know what you're doing -# num_dilations=3 #num iterations of images erosion, don't mess with unless you know what you're doing -# process_args=np.array([num_erosions,num_dilations]) - -threshold = 1.5 -# num iterations of 3d image stack dilations, don't mess with unless you know what you're doing -ome_dilation_iter = 1 - -#thresholds for grains in reconstructions -# only use orientations from grains with completnesses ABOVE this threshold -comp_thresh = 0.0 -chi2_thresh = 1.0 # only use orientations from grains BELOW this chi^2 - - -######reconstruction parameters -ome_range_deg = [(0., 359.75)] # degrees - -use_mask = True -#Mask info, used if use_mask=True -mask_data_file = '/LOC/tomo_mask.npz' -# this is generally the difference in y motor positions between the tomo and nf layer (tomo_motor_z-nf_motor_z), needed for registry -mask_vert_offset = -0.3 - -#these values will be used if no mask data is provided (use_mask=False) (no tomography) -# cross sectional to reconstruct (should be at least 20%-30% over sample width) -cross_sectional_dim = 1.35 -voxel_spacing = 0.005 # in mm, voxel spacing for the near field reconstruction - - -##vertical (y) reconstruction voxel bounds in mm, ALWAYS USED REGARDLESS OF MASK -# (mm) if bounds are equal, a single layer is produced -v_bnds = [-0.085, 0.085] - -beam_stop_y_cen = 0.0 # mm, measured from the origin of the detector paramters -beam_stop_width = 0.6 # mm, width of the beam stop verticallu -beam_stop_parms = np.array([beam_stop_y_cen, beam_stop_width]) - -##### multiprocessing controller parameters -check = None -limit = None -generate = None -ncpus = mp.cpu_count() -# chunksize for multiprocessing, don't mess with unless you know what you're doing -chunk_size = 500 - -RAM_set = False # if True, manually set max amount of ram -max_RAM = 256 # only used if RAM_set is true. in GB - -#### debug plotting -output_plot_check = True - - -#============================================================================== -# %% LOAD GRAIN AND EXPERIMENT DATA -#============================================================================== - -experiment, nf_to_ff_id_map = nfutil.gen_trial_exp_data(grain_out_file, det_file, - mat_file, mat_name, max_tth, - comp_thresh, chi2_thresh, misorientation_bnd, - misorientation_spacing, ome_range_deg, - nframes, beam_stop_parms) - - -#============================================================================== -# %% LOAD / GENERATE TEST DATA -#============================================================================== - - -if use_mask: - mask_data = np.load(mask_data_file) - - mask_full = mask_data['mask'] - Xs_mask = mask_data['Xs'] - Ys_mask = mask_data['Ys']-(mask_vert_offset) - Zs_mask = mask_data['Zs'] - voxel_spacing = mask_data['voxel_spacing'] - - # need to think about how to handle a single layer in this context - tomo_layer_centers = np.squeeze(Ys_mask[:, 0, 0]) - above = np.where(tomo_layer_centers >= v_bnds[0]) - below = np.where(tomo_layer_centers < v_bnds[1]) - - in_bnds = np.intersect1d(above, below) - - mask = mask_full[in_bnds] - Xs = Xs_mask[in_bnds] - Ys = Ys_mask[in_bnds] - Zs = Zs_mask[in_bnds] - - test_crds_full = np.vstack([Xs.flatten(), Ys.flatten(), Zs.flatten()]).T - n_crds = len(test_crds_full) - - to_use = np.where(mask.flatten())[0] - - -else: - test_crds_full, n_crds, Xs, Ys, Zs = nfutil.gen_nf_test_grid( - cross_sectional_dim, v_bnds, voxel_spacing) - to_use = np.arange(len(test_crds_full)) - -test_crds = test_crds_full[to_use, :] - - -#============================================================================== -# %% NEAR FIELD - MAKE MEDIAN DARK -#============================================================================== - -dark = nfutil.gen_nf_dark(data_folder, img_nums, num_for_dark, experiment.nrows, - experiment.ncols, dark_type='median', num_digits=6, stem=stem) - -#============================================================================== -# %% NEAR FIELD - LOAD IMAGE DATA AND PROCESS -#============================================================================== - -image_stack = nfutil.gen_nf_cleaned_image_stack(data_folder, img_nums, dark, experiment.nrows, experiment.ncols, process_type=process_type, - process_args=process_args, threshold=img_threshold, ome_dilation_iter=ome_dilation_iter, num_digits=6, stem=stem) - - -#============================================================================== -# %% NEAR FIELD - splitting -#============================================================================== - -if RAM_set is True: - RAM = max_RAM * 1e9 -else: - RAM = psutil.virtual_memory().available # in GB -# # turn into number of bytes - -RAM_to_use = 0.75 * RAM - -n_oris = len(nf_to_ff_id_map) -n_voxels = len(test_crds) - -bits_for_arrays = 64*n_oris*n_voxels + 192 * \ - n_voxels # bits raw conf + bits voxel positions -bytes_for_array = bits_for_arrays/8. - -n_groups = np.floor(bytes_for_array/RAM_to_use) # number of full groups -leftover_voxels = np.mod(n_voxels, n_groups) - -print('Splitting data into %d groups with %d leftover voxels' % - (int(n_groups), int(leftover_voxels))) - - -grouped_voxels = n_voxels - leftover_voxels - -if n_groups != 0: - voxels_per_group = grouped_voxels/n_groups -#============================================================================== -# %% BUILD MP CONTROLLER -#============================================================================== - -# assume that if os has fork, it will be used by multiprocessing. -# note that on python > 3.4 we could use multiprocessing get_start_method and -# set_start_method for a cleaner implementation of this functionality. -multiprocessing_start_method = 'fork' if hasattr(os, 'fork') else 'spawn' - - -controller = nfutil.build_controller( - ncpus=ncpus, chunk_size=chunk_size, check=check, generate=generate, limit=limit) - - -#============================================================================== -# %% TEST ORIENTATIONS -#============================================================================== - -## Would be nice to pack the images, quant_and_clip will need to be edited if this is added, dcp 6.2.2021 -#packed_image_stack = nfutil.dilate_image_stack(image_stack, experiment,controller) - -print('Testing Orientations...') -#%% Test orientations in groups - -if n_groups == 0: - raw_confidence = nfutil.test_orientations( - image_stack, experiment, test_crds, controller, multiprocessing_start_method) - - del controller - - raw_confidence_full = np.zeros( - [len(experiment.exp_maps), len(test_crds_full)]) - - for ii in np.arange(raw_confidence_full.shape[0]): - raw_confidence_full[ii, to_use] = raw_confidence[ii, :] - -else: - grain_map_list = np.zeros(n_voxels) - confidence_map_list = np.zeros(n_voxels) - - # test voxels in groups - for group in range(int(n_groups)): - voxels_to_test = test_crds[int( - group) * int(voxels_per_group):int(group + 1) * int(voxels_per_group), :] - print('Calculating group %d' % group) - raw_confidence = nfutil.test_orientations( - image_stack, experiment, voxels_to_test, controller, multiprocessing_start_method) - print('Calculated raw confidence group %d' % group) - grain_map_group_list, confidence_map_group_list = nfutil.process_raw_confidence( - raw_confidence, id_remap=nf_to_ff_id_map, min_thresh=0.0) - - grain_map_list[int( - group) * int(voxels_per_group):int(group + 1) * int(voxels_per_group)] = grain_map_group_list - - confidence_map_list[int( - group) * int(voxels_per_group):int(group + 1) * int(voxels_per_group)] = confidence_map_group_list - del raw_confidence - - if leftover_voxels > 0: - #now for the leftover voxels - voxels_to_test = test_crds[int( - n_groups) * int(voxels_per_group):, :] - raw_confidence = nfutil.test_orientations( - image_stack, experiment, voxels_to_test, controller, multiprocessing_start_method) - grain_map_group_list, confidence_map_group_list = nfutil.process_raw_confidence( - raw_confidence, id_remap=nf_to_ff_id_map, min_thresh=0.0) - - grain_map_list[int( - n_groups) * int(voxels_per_group):] = grain_map_group_list - - confidence_map_list[int( - n_groups) * int(voxels_per_group):] = confidence_map_group_list - - #fix so that chunking will work with tomography - grain_map_list_full = np.zeros(Xs.shape[0]*Xs.shape[1]*Xs.shape[2]) - confidence_map_list_full = np.zeros(Xs.shape[0]*Xs.shape[1]*Xs.shape[2]) - - for jj in np.arange(len(to_use)): - grain_map_list_full[to_use[jj]] = grain_map_list[jj] - confidence_map_list_full[to_use[jj]] = confidence_map_list[jj] - - #reshape them - grain_map = grain_map_list_full.reshape(Xs.shape) - confidence_map = confidence_map_list_full.reshape(Xs.shape) - -del controller - -#============================================================================== -# %% POST PROCESS W WHEN TOMOGRAPHY HAS BEEN USED -#============================================================================== - -#note all masking is already handled by not evaluating specific points -grain_map, confidence_map = nfutil.process_raw_confidence( - raw_confidence_full, Xs.shape, id_remap=nf_to_ff_id_map, min_thresh=0.0) - -#============================================================================== -# %% SAVE PROCESSED GRAIN MAP DATA -#============================================================================== - -nfutil.save_nf_data(output_dir, output_stem, grain_map, confidence_map, - Xs, Ys, Zs, experiment.exp_maps, id_remap=nf_to_ff_id_map) - -#============================================================================== -# %% CHECKING OUTPUT -#============================================================================== - - -mat = Material(name=mat_name, material_file=mat_file, dmin=valWUnit( - 'lp', 'length', 0.05, 'nm'), kev=valWUnit('kev', 'energy', beam_energy, 'keV')) - -if matplot: - if output_plot_check: - nfutil.plot_ori_map(grain_map, confidence_map, - experiment.exp_maps, 0, mat, id_remap=nf_to_ff_id_map) diff --git a/scripts/tiffs_from_h5.py b/scripts/tiffs_from_h5.py deleted file mode 100644 index b1c326851..000000000 --- a/scripts/tiffs_from_h5.py +++ /dev/null @@ -1,64 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Wed Jul 31 18:30:53 2019 - -@author: bernier2 -""" - -import os - -import numpy as np - -from hexrd import imageseries - -from skimage import io - - -# dirs -working_dir = '/Users/Shared/APS/PUP_AFRL_Feb19' -image_dir = os.path.join(working_dir, 'image_data') - -samp_name = 'ceria_cal' -scan_number = 0 - -tif_file_template = samp_name + '_%06d-%s.tif' - -raw_data_dir_template = os.path.join( - image_dir, - 'raw_images_%s_%06d-%s.yml' -) -yml_string = """ -image-files: - directory: %s - files: "%s" - -options: - empty-frames: 0 - max-frames: 0 -meta: - panel: %s -""" - - -ims = imageseries.open( - os.path.join(image_dir, 'ceria_cal.h5'), - 'hdf5', - path='/imageseries' -) - -metadata = ims.metadata -det_keys = np.array(metadata['panels'], dtype=str) - -for i, det_key in enumerate(det_keys): - yml_file = open( - raw_data_dir_template % (samp_name, scan_number, det_key), - 'w' - ) - tiff_fname = tif_file_template % (scan_number, det_key) - print(yml_string % (image_dir, tiff_fname, det_key), - file=yml_file) - io.imsave( - os.path.join(image_dir, tiff_fname), - ims[i] - )