Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add option to set eta_step size for eta omega maps #573

Merged
merged 5 commits into from
Dec 4, 2023
Merged

Add option to set eta_step size for eta omega maps #573

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
cd866e7
Add option to set eta_step size for eta omega maps
psavery Oct 27, 2023
4743958
Fix branch cut logic for smaller eta steps
psavery Nov 13, 2023
5a65d87
Remove histogramming in _generate_ring_params()
psavery Nov 29, 2023
a7c79e5
Fix logic for eta step
Dec 4, 2023
5da6bd9
Extract logic for fast_histogram into a function
psavery Dec 4, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions hexrd/config/findorientations.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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',
Expand Down
1 change: 1 addition & 0 deletions hexrd/findorientations.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
96 changes: 24 additions & 72 deletions hexrd/instrument/hedm_instrument.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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',
Expand Down Expand Up @@ -2666,68 +2666,26 @@ 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]
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]
)
pixel_etas = peta[pixel_ids]
reta_hist = histogram(pixel_etas, eta_edges)
bins_on_detector = np.where(reta_hist)[0]

# ring arc lenght on panel
arc_length = angularDifference(
eta_edges[reta_bin_idx[0]],
eta_edges[reta_bin_idx[-1]]
)
return pixel_etas, eta_edges, pixel_ids, bins_on_detector

# 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]
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
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):
Expand All @@ -2748,17 +2706,11 @@ def _run_histograms(rows, ims, tth_ranges, ring_maps, ring_params, threshold):
continue

# Unpack the params
retas, eta_bins, rtth_idx, reta_idx = params

if fast_histogram:
result = histogram1d(retas, len(eta_bins) - 1,
(eta_bins[0], eta_bins[-1]),
weights=image[rtth_idx])
else:
result, _ = histogram1d(retas, bins=eta_bins,
weights=image[rtth_idx])
pixel_etas, eta_edges, pixel_ids, bins_on_detector = params
result = histogram(pixel_etas, eta_edges, weights=image[pixel_ids])

this_map[i_row, reta_idx] = 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,
Expand Down