Many-grain density support (dense/streaky frames)

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:

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

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:

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,

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,

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",

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