diff --git a/hexrd/calibration/__init__.py b/hexrd/calibration/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/hexrd/calibration/instrument_calibrator.py b/hexrd/calibration/instrument_calibrator.py deleted file mode 100644 index 485b85141..000000000 --- a/hexrd/calibration/instrument_calibrator.py +++ /dev/null @@ -1,72 +0,0 @@ -import numpy as np -from scipy.optimize import leastsq, least_squares - - -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, - conv_tol=1e-4, fit_tth_tol=None, - 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 - - delta_r = np.inf - step_successful = True - while delta_r > conv_tol and step_successful: - data_dict = calib_class._extract_powder_lines( - fit_tth_tol=fit_tth_tol) - - # 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)) - else: - print('no improvement in residual!!!') - step_successful = False - return x1 diff --git a/hexrd/calibration/laue_calibrator.py b/hexrd/calibration/laue_calibrator.py deleted file mode 100755 index b1eb78da4..000000000 --- a/hexrd/calibration/laue_calibrator.py +++ /dev/null @@ -1,376 +0,0 @@ -import numpy as np - -from hexrd.matrixutil import findDuplicateVectors -from hexrd.fitting import fitpeak - - -class LaueCalibrator(object): - def __init__(self, instr, plane_data, img_dict, flags, crystal_params, - min_energy=5., max_energy=25., tth_tol=2.0, eta_tol=2.0, - 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.makeNew() # need a copy to munge - self._plane_data.wavelength = max_energy # force - self._plane_data.exclusions = None # don't need exclusions for Laue - self._img_dict = img_dict - self._params = np.asarray(crystal_params, dtype=float) - self._full_params = np.hstack( - [self._instr.calibration_parameters, self._params] - ) - assert len(flags) == len(self._full_params), \ - "flags must have %d elements" % len(self._full_params) - self._flags = flags - - # !!! scalar cutoffs for now - # TODO: make a list or spectral input - self._min_energy = min_energy - self._max_energy = max_energy - - # for polar interpolation - self._tth_tol = tth_tol - self._eta_tol = eta_tol - - # for peak fitting - # ??? fitting only, or do alternative peak detection? - self._pktype = pktype - - @property - def npi(self): - return len(self._instr.calibration_parameters) - - @property - def instr(self): - return self._instr - - @property - def plane_data(self): - self._plane_data.wavelength = self._instr.beam_energy - return self._plane_data - - @property - def img_dict(self): - return self._img_dict - - @property - def min_energy(self): - return self._min_energy - - @min_energy.setter - def min_energy(self, x): - assert np.isscalar(x), "min_energy must be a scalar value" - self._min_energy = x - - @property - def max_energy(self): - return self._max_energy - - @max_energy.setter - def max_energy(self, x): - assert np.isscalar(x), "max_energy must be a scalar value" - self._max_energy = x - - @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 params(self): - return self._params - - @params.setter - def params(self, x): - x = np.atleast_1d(x) - if len(x) != len(self.plane_data.lparms): - raise RuntimeError("params must have %d elements" - % len(self.plane_data.lparms)) - self._params = x - self._plane_data.lparms = x - - @property - def full_params(self): - return self._full_params - - @property - def npe(self): - return len(self._params) - - @property - def flags(self): - return self._flags - - @flags.setter - def flags(self, x): - x = np.atleast_1d(x) - nparams_instr = len(self.instr.calibration_parameters) - nparams_extra = len(self.params) - nparams = nparams_instr + nparams_extra - if len(x) != nparams: - raise RuntimeError("flags must have %d elements" % nparams) - self._flags = np.asarrasy(x, dtype=bool) - self._instr.calibration_flags = self._flags[:nparams_instr] - - @property - def pktype(self): - return self._pktype - - @pktype.setter - def pktype(self, x): - """ - currently only 'gaussian', 'lorentzian, 'pvoigt' or 'split_pvoigt' - """ - assert x in ['gaussian', 'lorentzian', 'pvoigt', 'split_pvoigt'], \ - "pktype '%s' not understood" - self._pktype = x - - # ========================================================================= - # METHODS - # ========================================================================= - - def _extract_peaks(self): - refl_patches = dict.fromkeys(self.instr.detectors) - for ip_key, ip in self.instr.detectors.items(): - patches = xrdutil.make_reflection_patches( - ip.config_dict(chi=self.instr.chi, tvec=self.instr.tvec, - beam_vector=self.instr.beam_vector), - valid_angs[ip_key], ip.angularPixelSize(valid_xy[ip_key]), - rmat_c=xfcapi.makeRotMatOfExpMap(expmap_c), - tth_tol=tth_tol, eta_tol=eta_tol, - npdiv=npdiv, quiet=True) - refl_patches[ip_key] = list(patches) - - -"""for labeling""" -labelStructure = ndimage.generate_binary_structure(2, 1) -refl_dict = dict.fromkeys(instr.detectors) -for ip_key, ip in instr.detectors.items(): - reflInfoList = [] - img = self.img_dict[ip_key] - native_area = ip.pixel_area - num_patches = len(refl_patches[ip_key]) - meas_xy = np.nan*np.ones((num_patches, 2)) - meas_angs = np.nan*np.ones((num_patches, 2)) - for iRefl, patch in enumerate(refl_patches[ip_key]): - # check for overrun - irow = patch[-1][0] - jcol = patch[-1][1] - if np.any([irow < 0, irow >= ip.rows, jcol < 0, jcol >= ip.cols]): - continue - if not np.all(ip.clip_to_panel(np.vstack([patch[1][0].flatten(), - patch[1][1].flatten()]).T)[1] - ): - continue - # use nearest interpolation - spot_data = img[irow, jcol] * patch[3] * npdiv**2 / native_area - fsd = snip2d_py3.snip2d(spot_data, w=5, order=2, numiter=2) - spot_data -= fsd - - spot_data -= np.amin(spot_data) - patch_size = spot_data.shape - - sigmax = np.round(0.5*np.min(spot_data.shape)/sigma_to_FWHM) - - # optional gaussian smoothing - if do_smoothing: - spot_data = filters.gaussian(spot_data, stdev) - - if use_blob_detection: - spot_data_scl = 2.*spot_data/np.max(spot_data) - 1. - - # Compute radii in the 3rd column. - blobs_log = blob_log(spot_data_scl, - min_sigma=5, - max_sigma=max(sigmax, 6), - num_sigma=int(sigmax-5), - threshold=blob_thresh, overlap=0.1) - numPeaks = len(blobs_log) - else: - labels, numPeaks = ndimage.label( - spot_data > np.percentile(spot_data, 99), - structure=labelStructure - ) - slabels = np.arange(1, numPeaks + 1) - tth_edges = patch[0][0][0, :] - eta_edges = patch[0][1][:, 0] - delta_tth = tth_edges[1] - tth_edges[0] - delta_eta = eta_edges[1] - eta_edges[0] - if numPeaks > 0: - peakId = iRefl - if use_blob_detection: - coms = blobs_log[:, :2] - else: - coms = np.array( - ndimage.center_of_mass( - spot_data, labels=labels, index=slabels - ) - ) - if numPeaks > 1: - # - center = np.r_[spot_data.shape]*0.5 - com_diff = coms - np.tile(center, (numPeaks, 1)) - closest_peak_idx = np.argmin(np.sum(com_diff**2, axis=1)) - # - else: - closest_peak_idx = 0 - pass # end multipeak conditional - # - coms = coms[closest_peak_idx] - ''' - # - fig, ax = plt.subplots() - ax.imshow(spot_data, cmap=cm.hot, interpolation='none') - ax.plot(blobs_log[:, 1], blobs_log[:, 0], 'c+') - ax.plot(coms[1], coms[0], 'ms', mfc='none') - fig.suptitle("panel %s, reflection %d" % (ip_key, iRefl)) - print("%s, %d, (%.2f, %.2f), (%d, %d)" - % (ip_key, iRefl, coms[0], coms[1], - patch_size[0], patch_size[1])) - # - ''' - assert(coms.ndim == 1), "oops" - # - if fit_peaks: - sigm = 0.2*np.min(spot_data.shape) - if use_blob_detection: - sigm = min(blobs_log[closest_peak_idx, 2], sigm) - y0, x0 = coms.flatten() - ampl = float(spot_data[int(y0), int(x0)]) - # y0, x0 = 0.5*np.array(spot_data.shape) - # ampl = np.max(spot_data) - a_par = c_par = 0.5/float(sigm**2) - b_par = 0. - bgx = bgy = 0. - bkg = np.min(spot_data) - params = [ampl, a_par, b_par, c_par, x0, y0, bgx, bgy, bkg] - # - result = leastsq(gaussian_2d, params, args=(spot_data)) - # - fit_par = result[0] - # - coms = np.array([fit_par[5], fit_par[4]]) - print("%s, %d, (%.2f, %.2f), (%d, %d)" - % (ip_key, iRefl, coms[0], coms[1], - patch_size[0], patch_size[1])) - row_cen = 0.1*patch_size[0] - col_cen = 0.1*patch_size[1] - # - if np.any( - [coms[0] < row_cen, coms[0] >= patch_size[0] - row_cen, - coms[1] < col_cen, coms[1] >= patch_size[1] - col_cen] - ): - continue - if (fit_par[0] < 1.): - continue - # - if plot_fits: - fit_g = ( - gaussian_2d(fit_par, spot_data) + spot_data.flatten() - ).reshape(spot_data.shape) - fig, ax = plt.subplots(1, 3, sharex=True, sharey=True) - fig.suptitle("panel %s, reflection %d" % (ip_key, iRefl)) - ax[0].imshow(spot_data, cmap=cm.hot, interpolation='none') - ax[0].plot(fit_par[4], fit_par[5], 'c+') - ax[1].imshow(fit_g, cmap=cm.hot, interpolation='none') - ax[2].imshow(spot_data - fit_g, - cmap=cm.hot, interpolation='none') - - # intensities - spot_intensity, int_err = nquad( - gaussian_2d_int, - [[0., 2.*y0], [0., 2.*x0]], - args=fit_par) - pass - com_angs = np.hstack([ - tth_edges[0] + (0.5 + coms[1])*delta_tth, - eta_edges[0] + (0.5 + coms[0])*delta_eta - ]) - - # grab intensities - if not fit_peaks: - if use_blob_detection: - spot_intensity = 10 - max_intensity = 10 - else: - spot_intensity = np.sum( - spot_data[labels == slabels[closest_peak_idx]] - ) - max_intensity = np.max( - spot_data[labels == slabels[closest_peak_idx]] - ) - else: - max_intensity = np.max(spot_data) - # need xy coords - cmv = np.atleast_2d(np.hstack([com_angs, omega])) - gvec_c = xfcapi.angles_to_gvec( - cmv, - chi=instr.chi, - rmat_c=rmat_c, - beam_vec=instr.beam_vector) - new_xy = xfcapi.gvec_to_xy( - gvec_c, - ip.rmat, rmat_s, rmat_c, - ip.tvec, instr.tvec, tvec_c, - beam_vec=instr.beam_vector) - meas_xy[iRefl, :] = new_xy - meas_angs[iRefl, :] = com_angs - else: - peakId = -999 - # - spot_intensity = np.nan - max_intensity = np.nan - pass - - reflInfoList.append([peakId, valid_hkls[ip_key][:, iRefl], - (spot_intensity, max_intensity), - valid_energy[ip_key][iRefl], - valid_angs[ip_key][iRefl, :], - meas_angs[iRefl, :], - meas_xy[iRefl, :]]) - pass - reflInfo = np.array( - [tuple(i) for i in reflInfoList], - dtype=__reflInfo_dtype) - refl_dict[ip_key] = reflInfo - - - def _update_from_reduced(self, reduced_parameters): - # first update instrument from input parameters - full_params = np.asarray(self.full_params) - full_params[self.flags] = reduced_params - - self.instr.update_from_parameter_list(full_params[:self.npi]) - self.params = full_params[self.npi:] - - - def model(self, reduced_parameters): - - self._update_from_reduced(reduced_parameters) - - # Laue pattern simulation - laue_simulation = self.instr.simulate_laue_pattern( - self.plane_data, - minEnergy=self.min_energy, - maxEnergy=self.max_energy, - grain_params=grain_params - ) - - # parse simulation results - valid_xy, valid_hkls, valid_energy, valid_angs = \ - parse_laue_simulation(laue_simulation) - - return valid_xy diff --git a/hexrd/calibration/plot_calibration.py b/hexrd/calibration/plot_calibration.py deleted file mode 100644 index 19cb14777..000000000 --- a/hexrd/calibration/plot_calibration.py +++ /dev/null @@ -1,41 +0,0 @@ -import numpy as np -import matplotlib.gridspec as gridspec -from matplotlib import pyplot as plt - - -data_dict = pc._extract_powder_lines(fit_tth_tol=1.0) - - -# %% sample plot to check fit line poistions ahead of fitting -frows = int(np.ceil(np.sqrt(instr.num_panels))) -fcols = int(np.floor(np.sqrt(instr.num_panels))) -fig, ax = plt.subplots(frows, fcols) -fig_row, fig_col = np.unravel_index(np.arange(instr.num_panels), - (frows, fcols)) - -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, tth_ranges = 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/hexrd/calibration/powder_calibrator.py b/hexrd/calibration/powder_calibrator.py deleted file mode 100644 index e624b1e8a..000000000 --- a/hexrd/calibration/powder_calibrator.py +++ /dev/null @@ -1,272 +0,0 @@ -import numpy as np - -from hexrd.matrixutil import findDuplicateVectors -from hexrd.fitting import fitpeak - - -class PowderCalibrator(object): - def __init__(self, instr, plane_data, img_dict, flags, - 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.makeNew() # need a copy to munge - self._plane_data.wavelength = self._instr.beam_energy # force - self._img_dict = img_dict - self._params = np.asarray(self.plane_data.lparms, dtype=float) - self._full_params = np.hstack( - [self._instr.calibration_parameters, self._params] - ) - assert len(flags) == len(self._full_params), \ - "flags must have %d elements" % len(self._full_params) - self._flags = flags - - # 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 npi(self): - return len(self._instr.calibration_parameters) - - @property - def instr(self): - return self._instr - - @property - def plane_data(self): - self._plane_data.wavelength = self._instr.beam_energy - return self._plane_data - - @property - def 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 params(self): - return self._params - - @params.setter - def params(self, x): - x = np.atleast_1d(x) - if len(x) != len(self.plane_data.lparms): - raise RuntimeError("params must have %d elements" - % len(self.plane_data.lparms)) - self._params = x - self._plane_data.lparms = x - - @property - def full_params(self): - return self._full_params - - @property - def npe(self): - return len(self._params) - - @property - def flags(self): - return self._flags - - @flags.setter - def flags(self, x): - x = np.atleast_1d(x) - nparams_instr = len(self.instr.calibration_parameters) - nparams_extra = len(self.params) - nparams = nparams_instr + nparams_extra - if len(x) != nparams: - raise RuntimeError("flags must have %d elements" % nparams) - self._flags = np.asarrasy(x, dtype=bool) - self._instr.calibration_flags = self._flags[:nparams_instr] - - @property - def pktype(self): - return self._pktype - - @pktype.setter - def pktype(self, x): - """ - currently only 'gaussian', 'lorentzian, 'pvoigt' or 'split_pvoigt' - """ - assert x in ['gaussian', 'lorentzian', 'pvoigt', 'split_pvoigt'], \ - "pktype '%s' not understood" - 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, fit_tth_tol=None): - """ - 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! - """ - if fit_tth_tol is None: - fit_tth_tol = self.tth_tol/4. - fit_tth_tol = np.radians(fit_tth_tol) - # ideal tth - wlen = self.instr.beam_wavelength - dsp_ideal = np.atleast_1d(self.plane_data.getPlaneSpacings()) - hkls_ref = self.plane_data.hkls.T - dsp0 = [] - hkls = [] - for idx in self.plane_data.getMergedRanges()[0]: - if len(idx) > 1: - eqv, uidx = findDuplicateVectors(np.atleast_2d(dsp_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] - dsp0.append(dsp_ideal[uidx]) - hkls.append(hkls_ref[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) - - # peak profile fitting - 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] - tth_pred = 2.*np.arcsin(0.5*wlen/dsp0[i_ring]) - center_err = abs(tth_meas - tth_pred) - if p[0] < 0.01 or center_err > fit_tth_tol: - tmp.append(np.empty((0, 5))) - continue - xy_meas = panel.angles_to_cart([[tth_meas, eta_ref], ]) - - # distortion - if panel.distortion is not None: - xy_meas = panel.distortion.apply_inverse(xy_meas) - - # cat results - tmp.append( - np.hstack( - [xy_meas.squeeze(), - tth_meas, - hkls[i_ring], - dsp0[i_ring], - eta_ref] - ) - ) - pass - if len(tmp) == 0: - rhs[det_key].append(np.empty((0, 5))) - else: - rhs[det_key].append(np.vstack(tmp)) - pass - pass - return rhs - - def _evaluate(self, reduced_params, data_dict, output='residual'): - """ - """ - # need this for dsp - bmat = self.plane_data.latVecOps['B'] - - # first update instrument from input parameters - full_params = np.asarray(self.full_params) - full_params[self.flags] = reduced_params - - self.instr.update_from_parameter_list(full_params[:self.npi]) - self.params = full_params[self.npi:] - - wlen = self.instr.beam_wavelength - - # build residual - retval = [] - for det_key, panel in self.instr.detectors.items(): - pdata = np.vstack(data_dict[det_key]) - if len(pdata) > 0: - hkls = pdata[:, 3:6] - gvecs = np.dot(hkls, bmat.T) - dsp0 = 1./np.sqrt(np.sum(gvecs*gvecs, axis=1)) - # dsp0 = pdata[:, -2] - eta0 = pdata[:, -1] - - # derive reference tth - tth0 = 2.*np.arcsin(0.5*wlen/dsp0) - calc_xy = panel.angles_to_cart(np.vstack([tth0, eta0]).T) - - # distortion if applicable, from ideal --> warped - if panel.distortion is not None: - calc_xy = panel.distortion.apply_inverse(calc_xy) - - if output == 'residual': - retval.append( - (pdata[:, :2].flatten() - calc_xy.flatten()) - ) - elif output == 'model': - retval.append( - calc_xy.flatten() - ) - else: - raise RuntimeError("unrecognized output flag '%s'" - % output) - else: - continue - return np.hstack(retval) - - def residual(self, reduced_params, data_dict): - return self._evaluate(reduced_params, data_dict) - - def model(self, reduced_params, data_dict): - return self._evaluate(reduced_params, data_dict, output='model')