From 047e6bcf827bd1a2eae47bc4f0c1872f82c56337 Mon Sep 17 00:00:00 2001 From: Hemant Sharma Date: Thu, 23 Jul 2026 11:12:01 -0500 Subject: [PATCH 1/2] Many-grain density: robust blur sizing, GaussSigmaMax cap, 4M match ceiling MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Enables indexing of dense, streaky, many-grain frames (hundreds of overlapping grains per image) that the pipeline previously mishandled. Validated against multi-grain synthetic ground truth during the July 2026 34-ID-E beam time. Three coupled changes: 1. laue_index/preprocess.py — robust spot spacing. The auto blur sigma was sized from the ABSOLUTE minimum neighbour distance, which collapses on dense frames: fragments of a single streak sit a few px apart, shrinking sigma to ~1 px and leaving the 0.4-deg orientation grid unable to reach spots up to ~18 px from a grid seed (recovery dropped to ~15% on a 150-grain synthetic frame). Fragments of the same streak SHOULD merge — they are one grain's reflection — so size the blur from the 10th percentile of neighbour distances instead of the minimum. This changes the default preprocessed output; the golden fixtures are updated to match. 2. GaussSigmaMax parameter (laue_stream_utils.py + preprocess.py) — optional cap on the auto sigma. On dense frames a large blur lights so much of the detector that chance matches flood the coarse gate and bury true grains; iterative-peel drivers start tight and widen per pass. Default 0.0 = no cap. 3. MAX_MATCHES 100k -> 4M (LaueMatchingGPU.cu, LaueMatchingGPUStream.cu). The kernel appends matches in arrival order, so overflow discards candidates ARBITRARILY — true grains included — before the top-score merge can rank them. 100k overflowed on many-grain frames. 4M x 12 B x MAX_STREAMS is < 200 MB on device and pinned host. Also in laue_stream_utils.py: read_solutions tolerates torn/truncated rows (daemon terminated mid-write on a dense frame), dropping malformed rows with a warning rather than failing the whole parse. Char suite and the streaming regression tests pass with these changes. --- laue_index/preprocess.py | 15 ++++++++++++++- scripts/laue_stream_utils.py | 19 ++++++++++++++++++- src/LaueMatchingGPU.cu | 6 +++++- src/LaueMatchingGPUStream.cu | 6 +++++- tests/golden/config_lsu_dict.json | 1 + tests/golden/preprocess_real_frame.json | 4 ++-- 6 files changed, 45 insertions(+), 6 deletions(-) diff --git a/laue_index/preprocess.py b/laue_index/preprocess.py index 5c4e1e3..b7435f2 100644 --- a/laue_index/preprocess.py +++ b/laue_index/preprocess.py @@ -216,7 +216,14 @@ def calculate_gaussian_sigma( tree = cKDTree(coords) # query k=2: first neighbor is the point itself (distance 0), second is nearest dists, _ = tree.query(coords, k=2) - min_px_dist = float(dists[:, 1].min()) + # Robust spot spacing: the absolute minimum pair distance collapses on + # dense/streaky frames (fragments of ONE streak sit a few px apart), which + # shrank sigma to ~1 px and left the 0.4-deg orientation grid unable to + # reach spots up to ~18 px from a grid seed (recovery dropped to ~15% on a + # 150-grain synthetic frame). Fragments of the same streak SHOULD merge -- + # they are one grain's reflection -- so size the blur from the 10th + # percentile of neighbour distances instead of the minimum. + min_px_dist = float(np.percentile(dists[:, 1], 10)) if pixel_size > 0 and distance > 0: delta = (distance * np.tan(np.radians(orient_spacing))) / pixel_size @@ -335,6 +342,12 @@ def preprocess_image( distance=cfg["distance"], orient_spacing=cfg["orientation_spacing"], ) + # Optional cap (param GaussSigmaMax): on DENSE frames a large blur lights + # so much of the detector that chance matches flood the coarse gate and + # bury true grains; iterative-peel drivers start tight and widen per pass. + smax = float(cfg.get("gauss_sigma_max", 0.0) or 0.0) + if smax > 0: + sigma = min(sigma, smax) blurred = ndimg.gaussian_filter(filt_img.astype(np.float64), sigma) if return_intermediates: diff --git a/scripts/laue_stream_utils.py b/scripts/laue_stream_utils.py index 0dbdb98..5b1385a 100755 --- a/scripts/laue_stream_utils.py +++ b/scripts/laue_stream_utils.py @@ -123,6 +123,7 @@ def setup_logger( "max_angle": 2.0, "max_laue_spots": 400, "orientation_spacing": 0.4, + "gauss_sigma_max": 0.0, # 0 = no cap on the auto matching-blur sigma # Files "background_file": "", "orientation_file": "orientations.bin", @@ -194,6 +195,8 @@ def parse_config(config_file: str) -> Dict[str, Any]: cfg["max_laue_spots"] = int(rest[0]) elif key == "OrientationSpacing": cfg["orientation_spacing"] = float(rest[0]) + elif key == "GaussSigmaMax": + cfg["gauss_sigma_max"] = float(rest[0]) elif key == "ThresholdMethod": cfg["threshold_method"] = rest[0].lower() elif key == "Threshold": @@ -344,7 +347,21 @@ def read_solutions(path: str) -> Tuple[np.ndarray, str]: try: data = np.loadtxt(path, skiprows=1) except ValueError: - data = np.genfromtxt(path, skip_header=1) + # Tolerate torn/truncated rows (e.g. daemon terminated mid-write on a + # dense frame): keep only rows with the modal column count. + import warnings + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + data = np.genfromtxt(path, skip_header=1, invalid_raise=False) + if data.ndim == 2: + good = ~np.isnan(data).any(axis=1) + if (~good).sum(): + logger = logging.getLogger(__name__) + logger.warning( + f"read_solutions: dropped {(~good).sum()} malformed row(s) " + f"from {path}" + ) + data = data[good] if data.size == 0: data = np.empty((0, 31)) elif data.ndim == 1: diff --git a/src/LaueMatchingGPU.cu b/src/LaueMatchingGPU.cu index a0811a0..0dec59a 100644 --- a/src/LaueMatchingGPU.cu +++ b/src/LaueMatchingGPU.cu @@ -40,7 +40,11 @@ inline void gpuAssert(cudaError_t code, const char *file, int line, } } -#define MAX_MATCHES 100000 // max matched orientations per image +#define MAX_MATCHES 4000000 // max matched orientations per image. +// Sized for DENSE (many-grain, streaky) frames: the kernel appends in +// arrival order, so an overflow here discards candidates ARBITRARILY -- +// true grains included -- before the top-score merge can rank them. +// 4M entries x 12 B x MAX_STREAMS is < 200 MB on device and pinned host. __global__ void compare(size_t nrPxX, size_t nrPxY, size_t nOr, size_t nrMaxSpots, float minInt, size_t minSps, diff --git a/src/LaueMatchingGPUStream.cu b/src/LaueMatchingGPUStream.cu index c98c31a..f06e079 100644 --- a/src/LaueMatchingGPUStream.cu +++ b/src/LaueMatchingGPUStream.cu @@ -67,7 +67,11 @@ inline void gpuAssert(cudaError_t code, const char *file, int line, } // ── CUDA kernel (same as LaueMatchingGPU.cu) ─────────────────────────── -#define MAX_MATCHES 100000 // max matched orientations per image +#define MAX_MATCHES 4000000 // max matched orientations per image. +// Sized for DENSE (many-grain, streaky) frames: the kernel appends in +// arrival order, so an overflow here discards candidates ARBITRARILY -- +// true grains included -- before the top-score merge can rank them. +// 4M entries x 12 B x MAX_STREAMS is < 200 MB on device and pinned host. __global__ void compare(size_t nrPxX, size_t nOr, size_t nrMaxSpots, float minInt, size_t minSps, uint16_t *oA, float *im, diff --git a/tests/golden/config_lsu_dict.json b/tests/golden/config_lsu_dict.json index 58ffadd..9a6eba3 100644 --- a/tests/golden/config_lsu_dict.json +++ b/tests/golden/config_lsu_dict.json @@ -8,6 +8,7 @@ "elo": 5.0, "enhance_contrast": false, "filter_radius": 101, + "gauss_sigma_max": 0.0, "gaussian_factor": 0.25, "h5_location": "/entry/data/data", "hkl_file": "valid_hkls_Ni.csv", diff --git a/tests/golden/preprocess_real_frame.json b/tests/golden/preprocess_real_frame.json index 0607df7..8690756 100644 --- a/tests/golden/preprocess_real_frame.json +++ b/tests/golden/preprocess_real_frame.json @@ -9,8 +9,8 @@ "sum": 0.0 }, "blurred": { - "max": 19063.720969, - "nonzero": 58360, + "max": 15405.58383, + "nonzero": 64671, "shape": [ 2048, 2048 From f92c05b0d338db5eabf94d8aad041a832a9bbc6b Mon Sep 17 00:00:00 2001 From: Hemant Sharma Date: Thu, 23 Jul 2026 11:12:01 -0500 Subject: [PATCH 2/2] GenerateSimulation: generate multi-grain synthetic frames MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add calc_recip_array / orientations_to_recips and accept a list of orientations, so a single frame can superpose many grains' patterns. This is the synthetic ground truth used to validate the many-grain density work — known orientations, known count, to measure recovery. --- scripts/GenerateSimulation.py | 51 +++++++++++++++++++++++++++++++---- 1 file changed, 46 insertions(+), 5 deletions(-) diff --git a/scripts/GenerateSimulation.py b/scripts/GenerateSimulation.py index 0f73d36..bdbd9ce 100755 --- a/scripts/GenerateSimulation.py +++ b/scripts/GenerateSimulation.py @@ -157,6 +157,46 @@ def load_orientations(self): sys.exit(1) +def calc_recip_array(lat, sg_num): + """Reciprocal-lattice basis B (3x3) from lattice parameters. + + Port of calcRecipArray in src/LaueMatchingHeaders.h (and GenerateHKLs.py), + so q = OM @ B @ hkl matches the C matcher for ALL crystal systems. The + previous `OM * astar` shortcut is only valid for cubic lattices and put + simulated spots in the wrong places for e.g. hexagonal Ti-alpha. + For cubic, B == astar * I, so cubic output is unchanged. + """ + from math import cos, sin, sqrt, pi as _pi + a, b, c, alpha, beta, gamma = [float(x) for x in lat] + d2r = _pi / 180.0 + rhomb = 1 if sg_num in [146, 148, 155, 160, 161, 166, 167] else 0 + ca, cb, cg, sg = cos(alpha*d2r), cos(beta*d2r), cos(gamma*d2r), sin(gamma*d2r) + phi = sqrt(1.0 - ca*ca - cb*cb - cg*cg + 2*ca*cb*cg) + pv = (2*_pi) / (a*b*c*phi) + zo = lambda v: 0.0 if abs(v) < 1e-11 else v + if rhomb == 0: + a0, a1, a2 = a, 0.0, 0.0 + b0, b1, b2 = b*cg, b*sg, 0.0 + c0, c1, c2 = c*cb, c*(ca-cb*cg)/sg, c*phi/sg + else: + p = sqrt(1.0 + 2*ca); q = sqrt(1.0 - ca) + pmq = (a/3.0)*(p-q); p2q = (a/3.0)*(p+2*q) + a0, a1, a2 = p2q, pmq, pmq + b0, b1, b2 = pmq, p2q, pmq + c0, c1, c2 = pmq, pmq, p2q + B = np.zeros((3, 3)) + B[0][0] = zo((b1*c2-b2*c1)*pv); B[1][0] = zo((b2*c0-b0*c2)*pv); B[2][0] = zo((b0*c1-b1*c0)*pv) + B[0][1] = zo((c1*a2-c2*a1)*pv); B[1][1] = zo((c2*a0-c0*a2)*pv); B[2][1] = zo((c0*a1-c1*a0)*pv) + B[0][2] = zo((a1*b2-a2*b1)*pv); B[1][2] = zo((a2*b0-a0*b2)*pv); B[2][2] = zo((a0*b1-a1*b0)*pv) + return B + + +def orientations_to_recips(orientations, params): + """Rows of flattened OMs -> rows of flattened OM @ B (crystal-correct).""" + B = calc_recip_array(params['latC'].split(), params['sgNum']) + return np.array([(om.reshape(3, 3) @ B).flatten() for om in orientations]) + + class DiffractionSimulator: """Simulate diffraction patterns based on orientations and configuration.""" @@ -345,9 +385,10 @@ def simulate(self, orientations): Returns: tuple: (image, position array) """ - # Scale orientations by astar - recips = orientations * self.params['astar'] - + # Rotate the crystal-correct reciprocal basis (NOT the cubic-only + # OM*astar shortcut; see calc_recip_array). + recips = orientations_to_recips(orientations, self.params) + # Process each orientation for grain_nr, recip in enumerate(recips): recip = recip.reshape((3, 3)) @@ -380,8 +421,8 @@ def save_results(self, output_file, orientations, provenance_inputs=None): with h5py.File(output_file, 'w') as hf: hf.create_dataset('/entry1/data/data', data=self.img) - # Reshape recips to store in file - recips = orientations * self.params['astar'] + # Reshape recips to store in file (crystal-correct basis) + recips = orientations_to_recips(orientations, self.params) hf.create_dataset('/entry1/recips', data=recips) # Store spot positions