|
| 1 | +import copy |
| 2 | +from pathlib import Path |
| 3 | + |
| 4 | +import h5py |
| 5 | +import numpy as np |
| 6 | + |
| 7 | +import pytest |
| 8 | + |
| 9 | +from hexrd.instrument.hedm_instrument import HEDMInstrument |
| 10 | +from hexrd.material.material import load_materials_hdf5, Material |
| 11 | +from hexrd.projections.polar import PolarView |
| 12 | + |
| 13 | + |
| 14 | +@pytest.fixture |
| 15 | +def eiger_examples_path(example_repo_path: Path) -> Path: |
| 16 | + return Path(example_repo_path) / 'eiger' |
| 17 | + |
| 18 | + |
| 19 | +@pytest.fixture |
| 20 | +def ceria_examples_path(eiger_examples_path: Path) -> Path: |
| 21 | + return eiger_examples_path / 'first_ceria' |
| 22 | + |
| 23 | + |
| 24 | +@pytest.fixture |
| 25 | +def eiger_instrument(ceria_examples_path: Path) -> HEDMInstrument: |
| 26 | + instr_path = ( |
| 27 | + ceria_examples_path / 'eiger_ceria_calibrated_composite.hexrd' |
| 28 | + ) |
| 29 | + with h5py.File(instr_path, 'r') as rf: |
| 30 | + return HEDMInstrument(rf) |
| 31 | + |
| 32 | + |
| 33 | +@pytest.fixture |
| 34 | +def ceria_material(ceria_examples_path: Path) -> Material: |
| 35 | + path = ceria_examples_path / 'ceria.h5' |
| 36 | + materials = load_materials_hdf5(path) |
| 37 | + return materials['CeO2'] |
| 38 | + |
| 39 | + |
| 40 | +@pytest.fixture |
| 41 | +def expected_simulated_powder_lines_results( |
| 42 | + test_data_dir: Path, |
| 43 | +) -> dict[str, dict]: |
| 44 | + path = test_data_dir / 'expected_simulated_powder_lines_results.npy' |
| 45 | + return np.load(path, allow_pickle=True).item() |
| 46 | + |
| 47 | + |
| 48 | +def test_simulate_powder_lines( |
| 49 | + eiger_instrument: HEDMInstrument, |
| 50 | + ceria_material: Material, |
| 51 | + expected_simulated_powder_lines_results: dict[str, dict], |
| 52 | +): |
| 53 | + instr = eiger_instrument |
| 54 | + material = copy.deepcopy(ceria_material) |
| 55 | + ref = expected_simulated_powder_lines_results |
| 56 | + |
| 57 | + pd = material.planeData |
| 58 | + |
| 59 | + pd.exclusions = None |
| 60 | + |
| 61 | + pd.tThWidth = np.radians(0.25) |
| 62 | + |
| 63 | + hkls = pd.getHKLs() |
| 64 | + tth = pd.getTTh() |
| 65 | + |
| 66 | + # There should be exactly 22 two HKLs generated |
| 67 | + assert len(tth) == 22 |
| 68 | + |
| 69 | + def hkl_idx(hkl): |
| 70 | + hkl = list(hkl) |
| 71 | + for i, hkl_ref in enumerate(hkls.tolist()): |
| 72 | + if hkl_ref == hkl: |
| 73 | + return i |
| 74 | + |
| 75 | + return None |
| 76 | + |
| 77 | + # Verify a few manual cases |
| 78 | + # The four corners should all have these HKLs |
| 79 | + corner_hkls = { |
| 80 | + (4, 2, 0): 9.583035640913781, |
| 81 | + # 333 and 511 are nearly identical. |
| 82 | + (3, 3, 3): 11.139041445220014, |
| 83 | + (5, 1, 1): 11.139041445220013, |
| 84 | + (6, 2, 0): 13.568350297337751, |
| 85 | + } |
| 86 | + corner_hkls = {k: np.radians(v) for k, v in corner_hkls.items()} |
| 87 | + |
| 88 | + # The four inner subpanels should *only* have these HKLs |
| 89 | + inner_hkls = { |
| 90 | + (1, 1, 1): 3.7078160949754677, |
| 91 | + (2, 0, 0): 4.281666411410814, |
| 92 | + } |
| 93 | + inner_hkls = {k: np.radians(v) for k, v in inner_hkls.items()} |
| 94 | + |
| 95 | + inner_dets = ['eiger_3_1', 'eiger_3_2', 'eiger_4_1', 'eiger_4_2'] |
| 96 | + corner_dets = ['eiger_0_0', 'eiger_0_3', 'eiger_7_0', 'eiger_7_3'] |
| 97 | + |
| 98 | + full_results = {} |
| 99 | + for det_key, panel in instr.detectors.items(): |
| 100 | + angs, xys, ranges = panel.make_powder_rings( |
| 101 | + pd, |
| 102 | + # Have to make a smaller delta eta to get all intersections, |
| 103 | + # since the Eiger subpanels are small. |
| 104 | + delta_eta=0.25, |
| 105 | + ) |
| 106 | + full_results[det_key] = { |
| 107 | + 'angs': angs, |
| 108 | + 'xys': xys, |
| 109 | + 'ranges': ranges, |
| 110 | + } |
| 111 | + |
| 112 | + valid_angs = [x for x in angs if x.size > 0] |
| 113 | + if det_key in inner_dets: |
| 114 | + for hkl, value in inner_hkls.items(): |
| 115 | + idx = hkl_idx(hkl) |
| 116 | + assert np.allclose(angs[idx][0][0], value) |
| 117 | + |
| 118 | + # Also verify that there are no other HKLs |
| 119 | + assert len(valid_angs) == len(inner_hkls) |
| 120 | + |
| 121 | + if det_key in corner_dets: |
| 122 | + for hkl, value in corner_hkls.items(): |
| 123 | + idx = hkl_idx(hkl) |
| 124 | + assert np.allclose(angs[idx][0][0], value) |
| 125 | + |
| 126 | + # There should be at least 9 valid HKLs on every corner detector |
| 127 | + assert len(valid_angs) >= 9 |
| 128 | + |
| 129 | + # Verify that a few of the HKLs were combined |
| 130 | + indices, ranges = pd.getMergedRanges() |
| 131 | + |
| 132 | + num_merged_hkls = sum([len(x) > 1 for x in indices]) |
| 133 | + assert num_merged_hkls == 5 |
| 134 | + |
| 135 | + # Verify that 3, 3, 3 and 5, 1, 1 were merged |
| 136 | + idx1 = hkl_idx([5, 1, 1]) |
| 137 | + idx2 = hkl_idx([3, 3, 3]) |
| 138 | + assert [idx1, idx2] in indices |
| 139 | + |
| 140 | + # Now just verify that everything matches what it was before... |
| 141 | + for det_key, det_results in full_results.items(): |
| 142 | + for entry_key, results in det_results.items(): |
| 143 | + ref_results = ref[det_key][entry_key] |
| 144 | + for i in range(len(results)): |
| 145 | + assert np.allclose(results[i], ref_results[i], equal_nan=True) |
| 146 | + |
| 147 | + |
| 148 | +def test_simulate_powder_pattern_image( |
| 149 | + eiger_instrument: HEDMInstrument, |
| 150 | + ceria_material: Material, |
| 151 | +): |
| 152 | + instr = eiger_instrument |
| 153 | + material = ceria_material |
| 154 | + |
| 155 | + img_dict = instr.simulate_powder_pattern( |
| 156 | + [material], |
| 157 | + noise='poisson', |
| 158 | + ) |
| 159 | + |
| 160 | + # Now do a warp to the polar view, create the lineout, and verify there |
| 161 | + # is intensity in a few places we'd expect. |
| 162 | + tth_range = [0.1, 14] |
| 163 | + eta_min = -180 |
| 164 | + eta_max = 180 |
| 165 | + pixel_size = (0.1, 0.1) |
| 166 | + pv = PolarView(tth_range, instr, eta_min, eta_max, pixel_size) |
| 167 | + img = pv.warp_image(img_dict, pad_with_nans=True, |
| 168 | + do_interpolation=True) |
| 169 | + |
| 170 | + lineout = img.mean(axis=0).filled(np.nan) |
| 171 | + |
| 172 | + # Now verify there is intensity at the predicted two theta values |
| 173 | + # Sort by structure factor and check the 4 most intense lines |
| 174 | + sf_sorting = np.argsort(-material.planeData.structFact) |
| 175 | + for i in range(4): |
| 176 | + tth_idx = sf_sorting[i] |
| 177 | + tth_val = np.degrees(material.planeData.getTTh()[tth_idx]) |
| 178 | + |
| 179 | + # Compute the index for this tth value in the lineout |
| 180 | + idx = int(np.floor((tth_val - tth_range[0]) / pixel_size[1])) |
| 181 | + |
| 182 | + # Verify we are at a maximum, and that the value is higher |
| 183 | + # than the background. |
| 184 | + assert lineout[idx] > lineout[0] |
| 185 | + assert lineout[idx] > lineout[idx - 1] |
| 186 | + assert lineout[idx] > lineout[idx + 1] |
0 commit comments