From cd866e7410174857ce561e6108ad5d06d249a673 Mon Sep 17 00:00:00 2001 From: Patrick Avery Date: Fri, 27 Oct 2023 12:35:06 -0500 Subject: [PATCH 1/5] Add option to set eta_step size for eta omega maps It was previously fixed to 0.25 degrees. However, the users may find it beneficial to change this. Smaller eta_step sizes might lead to less overlap between neighboring peaks in the eta omega maps, for example. Signed-off-by: Patrick Avery --- hexrd/config/findorientations.py | 6 ++++++ hexrd/findorientations.py | 1 + 2 files changed, 7 insertions(+) diff --git a/hexrd/config/findorientations.py b/hexrd/config/findorientations.py index 5326796b7..c66034bc8 100644 --- a/hexrd/config/findorientations.py +++ b/hexrd/config/findorientations.py @@ -223,6 +223,12 @@ def bin_frames(self): 'find_orientations:orientation_maps:bin_frames', default=1 ) + @property + def eta_step(self): + return self._cfg.get( + 'find_orientations:orientation_maps:eta_step', default=0.25 + ) + @property def file(self): temp = self._cfg.get('find_orientations:orientation_maps:file', diff --git a/hexrd/findorientations.py b/hexrd/findorientations.py index 76721972e..36f08f1cc 100755 --- a/hexrd/findorientations.py +++ b/hexrd/findorientations.py @@ -567,6 +567,7 @@ def generate_eta_ome_maps(cfg, hkls=None, save=True): eta_ome = instrument.GenerateEtaOmeMaps( imsd, cfg.instrument.hedm, plane_data, active_hkls=active_hklIDs, + eta_step=cfg.find_orientations.orientation_maps.eta_step, threshold=cfg.find_orientations.orientation_maps.threshold, ome_period=ome_period) From 4743958fae2e3f11796110b8dd331675d297ca19 Mon Sep 17 00:00:00 2001 From: Patrick Avery Date: Mon, 13 Nov 2023 16:37:44 -0600 Subject: [PATCH 2/5] Fix branch cut logic for smaller eta steps Smaller values in eta step were breaking the branch cut logic here. The problem is that the arc length was not always less than 1e-4 when the eta_step was small (such as 0.05). It was less than 1e-2, so still close to zero, but not small enough. Then on top of that, some of the ring break logic was not working. To fix both of those issues, we loosen the logic to instead check to see if nearly the whole two pi range is covered, and if there is a big gap. If there is, then we assume it is a branch cut. And we locate the break via the big gap. This worked properly for eta step values from 0.01 to 1 Signed-off-by: Patrick Avery --- hexrd/instrument/hedm_instrument.py | 39 +++++++++++------------------ 1 file changed, 15 insertions(+), 24 deletions(-) diff --git a/hexrd/instrument/hedm_instrument.py b/hexrd/instrument/hedm_instrument.py index 889591549..6c0d0545a 100644 --- a/hexrd/instrument/hedm_instrument.py +++ b/hexrd/instrument/hedm_instrument.py @@ -2688,33 +2688,24 @@ def _generate_ring_params(tthr, ptth, peta, eta_edges, delta_eta): [reta_idx, reta_idx[-1] + 1] ) + eta_bins = eta_edges[reta_bin_idx] - # ring arc lenght on panel - arc_length = angularDifference( - eta_edges[reta_bin_idx[0]], - eta_edges[reta_bin_idx[-1]] - ) + arc_length = eta_edges[reta_bin_idx][-1] - eta_edges[reta_bin_idx][0] + etas_span_2pi = (arc_length / (2 * np.pi)) > 0.99 - # Munge eta bins - # !!! need to work with the subset to preserve - # NaN values at panel extents! - # - # !!! MUST RE-MAP IF BRANCH CUT IS IN RANGE - # - # The logic below assumes that eta_edges span 2*pi to - # single precision - eta_bins = eta_edges[reta_bin_idx] - if arc_length < 1e-4: - # have branch cut in here - ring_gap = np.where( - reta_idx - - np.arange(len(reta_idx)) - )[0] - if len(ring_gap) > 0: - # have incomplete ring - eta_stop_idx = ring_gap[0] - eta_stop = eta_edges[eta_stop_idx] + if etas_span_2pi: + eta_bins_diff = np.diff(eta_bins) + has_big_gap = eta_bins_diff.max() > np.median(eta_bins_diff) * 2 + + if has_big_gap: + # If there is a big gap and the eta values span nearly the + # whole two pi range, assume we have a branch cut. We must remap. + + # Find the biggest gap. + eta_stop_idx = np.argmax(eta_bins_diff) + 1 + eta_stop = eta_bins[eta_stop_idx] new_period = np.cumsum([eta_stop, 2*np.pi]) + # remap retas = mapAngle(retas, new_period) tmp_bins = mapAngle( From 5a65d878a1ec4f2fd825f03d398484c3fc152afd Mon Sep 17 00:00:00 2001 From: Patrick Avery Date: Wed, 29 Nov 2023 12:35:19 -0600 Subject: [PATCH 3/5] Remove histogramming in _generate_ring_params() This change was suggested by @donald-e-boyce He mentioned that the whole histogramming should be unnecessary, and we should just simplify the code to not use it. This appears to work properly for the different `eta_step` sizes as before. Thanks Don! Signed-off-by: Patrick Avery --- hexrd/instrument/hedm_instrument.py | 56 ++++------------------------- 1 file changed, 6 insertions(+), 50 deletions(-) diff --git a/hexrd/instrument/hedm_instrument.py b/hexrd/instrument/hedm_instrument.py index 6c0d0545a..b839ec146 100644 --- a/hexrd/instrument/hedm_instrument.py +++ b/hexrd/instrument/hedm_instrument.py @@ -61,7 +61,6 @@ from hexrd import matrixutil as mutil from hexrd.transforms.xfcapi import ( angles_to_gvec, - angularDifference, gvec_to_xy, make_sample_rmat, makeRotMatOfExpMap, @@ -2675,50 +2674,7 @@ def _generate_ring_params(tthr, ptth, peta, eta_edges, delta_eta): # detect a branch cut. The histogram idx var # is the left-hand edges... retas = peta[rtth_idx] - if fast_histogram: - reta_hist = histogram1d( - retas, - len(eta_edges) - 1, - (eta_edges[0], eta_edges[-1]) - ) - else: - reta_hist, _ = histogram1d(retas, bins=eta_edges) - reta_idx = np.where(reta_hist)[0] - reta_bin_idx = np.hstack( - [reta_idx, - reta_idx[-1] + 1] - ) - eta_bins = eta_edges[reta_bin_idx] - - arc_length = eta_edges[reta_bin_idx][-1] - eta_edges[reta_bin_idx][0] - etas_span_2pi = (arc_length / (2 * np.pi)) > 0.99 - - if etas_span_2pi: - eta_bins_diff = np.diff(eta_bins) - has_big_gap = eta_bins_diff.max() > np.median(eta_bins_diff) * 2 - - if has_big_gap: - # If there is a big gap and the eta values span nearly the - # whole two pi range, assume we have a branch cut. We must remap. - - # Find the biggest gap. - eta_stop_idx = np.argmax(eta_bins_diff) + 1 - eta_stop = eta_bins[eta_stop_idx] - new_period = np.cumsum([eta_stop, 2*np.pi]) - - # remap - retas = mapAngle(retas, new_period) - tmp_bins = mapAngle( - eta_edges[reta_idx], new_period - ) - tmp_idx = np.argsort(tmp_bins) - reta_idx = reta_idx[np.argsort(tmp_bins)] - eta_bins = np.hstack( - [tmp_bins[tmp_idx], - tmp_bins[tmp_idx][-1] + delta_eta] - ) - - return retas, eta_bins, rtth_idx, reta_idx + return rtth_idx, retas, eta_edges def _run_histograms(rows, ims, tth_ranges, ring_maps, ring_params, threshold): @@ -2739,17 +2695,17 @@ def _run_histograms(rows, ims, tth_ranges, ring_maps, ring_params, threshold): continue # Unpack the params - retas, eta_bins, rtth_idx, reta_idx = params + rtth_idx, retas, eta_edges = params if fast_histogram: - result = histogram1d(retas, len(eta_bins) - 1, - (eta_bins[0], eta_bins[-1]), + result = histogram1d(retas, len(eta_edges) - 1, + (eta_edges[0], eta_edges[-1]), weights=image[rtth_idx]) else: - result, _ = histogram1d(retas, bins=eta_bins, + result, _ = histogram1d(retas, bins=eta_edges, weights=image[rtth_idx]) - this_map[i_row, reta_idx] = result + this_map[i_row, :] = result def _extract_detector_line_positions(iter_args, plane_data, tth_tol, From a7c79e5ad331c1db1284910e7a88ad4b403663f9 Mon Sep 17 00:00:00 2001 From: "donald e. boyce" Date: Mon, 4 Dec 2023 10:55:40 -0600 Subject: [PATCH 4/5] Fix logic for eta step --- hexrd/instrument/hedm_instrument.py | 37 +++++++++++++++++------------ 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/hexrd/instrument/hedm_instrument.py b/hexrd/instrument/hedm_instrument.py index b839ec146..c7054d237 100644 --- a/hexrd/instrument/hedm_instrument.py +++ b/hexrd/instrument/hedm_instrument.py @@ -61,6 +61,7 @@ from hexrd import matrixutil as mutil from hexrd.transforms.xfcapi import ( angles_to_gvec, + angularDifference, gvec_to_xy, make_sample_rmat, makeRotMatOfExpMap, @@ -2665,16 +2666,22 @@ def _generate_ring_params(tthr, ptth, peta, eta_edges, delta_eta): if not np.any(pixels_in_tthr): return None - # ???: faster to index with bool or use np.where, - # or recode in numba? - rtth_idx = np.where(pixels_in_tthr) + pixel_ids = np.where(pixels_in_tthr) # grab relevant eta coords using histogram - # !!!: This allows use to calculate arc length and - # detect a branch cut. The histogram idx var - # is the left-hand edges... - retas = peta[rtth_idx] - return rtth_idx, retas, eta_edges + pixel_etas = peta[pixel_ids] + if fast_histogram: + reta_hist = histogram1d( + pixel_etas, + len(eta_edges) - 1, + (eta_edges[0], eta_edges[-1]) + ) + else: + reta_hist, _ = histogram1d(pixel_etas, bins=eta_edges) + + bins_on_detector = np.where(reta_hist)[0] + + return pixel_etas, eta_edges, pixel_ids, bins_on_detector def _run_histograms(rows, ims, tth_ranges, ring_maps, ring_params, threshold): @@ -2695,17 +2702,17 @@ def _run_histograms(rows, ims, tth_ranges, ring_maps, ring_params, threshold): continue # Unpack the params - rtth_idx, retas, eta_edges = params - + pixel_etas, eta_edges, pixel_ids, bins_on_detector = params if fast_histogram: - result = histogram1d(retas, len(eta_edges) - 1, + result = histogram1d(pixel_etas, len(eta_edges) - 1, (eta_edges[0], eta_edges[-1]), - weights=image[rtth_idx]) + weights=image[pixel_ids]) else: - result, _ = histogram1d(retas, bins=eta_edges, - weights=image[rtth_idx]) + result, _ = histogram1d(pixel_etas, bins=eta_edges, + weights=image[pixel_ids]) - this_map[i_row, :] = result + # Note that this preserves nan values for bins not on the detector. + this_map[i_row, bins_on_detector] = result[bins_on_detector] def _extract_detector_line_positions(iter_args, plane_data, tth_tol, From 5da6bd95c7e486ec8e85c6077a969d351db02c59 Mon Sep 17 00:00:00 2001 From: Patrick Avery Date: Mon, 4 Dec 2023 11:11:31 -0600 Subject: [PATCH 5/5] Extract logic for fast_histogram into a function This is easier than the logic checks we were doing before. Signed-off-by: Patrick Avery --- hexrd/instrument/hedm_instrument.py | 40 ++++++++++++++--------------- 1 file changed, 19 insertions(+), 21 deletions(-) diff --git a/hexrd/instrument/hedm_instrument.py b/hexrd/instrument/hedm_instrument.py index c7054d237..c117f9903 100644 --- a/hexrd/instrument/hedm_instrument.py +++ b/hexrd/instrument/hedm_instrument.py @@ -61,7 +61,6 @@ from hexrd import matrixutil as mutil from hexrd.transforms.xfcapi import ( angles_to_gvec, - angularDifference, gvec_to_xy, make_sample_rmat, makeRotMatOfExpMap, @@ -649,8 +648,9 @@ def __init__(self, instrument_config=None, if det_buffer.ndim == 2: if det_buffer.shape != shape: msg = ( - f'Buffer shape for {det_id} ({det_buffer.shape}) ' - f'does not match detector shape ({shape})' + f'Buffer shape for {det_id} ' + f'({det_buffer.shape}) does not match ' + f'detector shape ({shape})' ) raise BufferShapeMismatchError(msg) else: @@ -1111,8 +1111,8 @@ def update_from_lmfit_parameter_list(self, params): for det_name, detector in self.detectors.items(): det = det_name.replace('-', '_') euler = np.r_[params[f'{det}_euler_z'].value, - params[f'{det}_euler_xp'].value, - params[f'{det}_euler_zpp'].value] + params[f'{det}_euler_xp'].value, + params[f'{det}_euler_zpp'].value] rmat = make_rmat_euler(np.radians(euler), 'zxz', @@ -2670,20 +2670,24 @@ def _generate_ring_params(tthr, ptth, peta, eta_edges, delta_eta): # grab relevant eta coords using histogram pixel_etas = peta[pixel_ids] - if fast_histogram: - reta_hist = histogram1d( - pixel_etas, - len(eta_edges) - 1, - (eta_edges[0], eta_edges[-1]) - ) - else: - reta_hist, _ = histogram1d(pixel_etas, bins=eta_edges) - + reta_hist = histogram(pixel_etas, eta_edges) bins_on_detector = np.where(reta_hist)[0] return pixel_etas, eta_edges, pixel_ids, bins_on_detector +def run_fast_histogram(x, bins, weights=None): + return histogram1d(x, len(bins) - 1, (bins[0], bins[-1]), + weights=weights) + + +def run_numpy_histogram(x, bins, weights=None): + return histogram1d(x, bins=bins, weights=weights)[0] + + +histogram = run_fast_histogram if fast_histogram else run_numpy_histogram + + def _run_histograms(rows, ims, tth_ranges, ring_maps, ring_params, threshold): for i_row in range(*rows): image = ims[i_row] @@ -2703,13 +2707,7 @@ def _run_histograms(rows, ims, tth_ranges, ring_maps, ring_params, threshold): # Unpack the params pixel_etas, eta_edges, pixel_ids, bins_on_detector = params - if fast_histogram: - result = histogram1d(pixel_etas, len(eta_edges) - 1, - (eta_edges[0], eta_edges[-1]), - weights=image[pixel_ids]) - else: - result, _ = histogram1d(pixel_etas, bins=eta_edges, - weights=image[pixel_ids]) + result = histogram(pixel_etas, eta_edges, weights=image[pixel_ids]) # Note that this preserves nan values for bins not on the detector. this_map[i_row, bins_on_detector] = result[bins_on_detector]