perl_d
78b2bb7a88
CI / lint (push) Skipped
CI / test (3.11) (push) Skipped
CI / test (3.12) (push) Skipped
CI / test (3.13) (push) Skipped
CI / lint (pull_request) Failing after 41s
CI / test (3.11) (pull_request) Skipped
CI / test (3.12) (pull_request) Skipped
CI / test (3.13) (pull_request) Skipped
695 lines
25 KiB
Python
695 lines
25 KiB
Python
from collections.abc import Callable
|
||
|
||
import numpy as np
|
||
from scipy import ndimage
|
||
|
||
from aarecommon.config.logger import setup_logger
|
||
from aarecommon.models.models import CrystalSize
|
||
from aarecommon.models.raster_grid import CenterOfMassModel, RasterGridRequest
|
||
|
||
logger = setup_logger("aareDAQ")
|
||
|
||
|
||
def identify_crystal_raster(result, r: RasterGridRequest) -> CenterOfMassModel | None:
|
||
images = result.images
|
||
if images and any(getattr(img, "spots_low_res", 0) for img in images):
|
||
logger.debug("Find image by maximum number of low resolution spots")
|
||
max_image = max(images, key=lambda img: img.spots_low_res)
|
||
logger.debug(f"Image with maximum spots_low_res: {max_image}")
|
||
logger.debug(f"Maximum spots_low_res value: {max_image.spots_low_res}")
|
||
logger.debug(f"Maximum image found at grid coordiantes {max_image.nx}, {max_image.ny}")
|
||
logger.debug(
|
||
f"Maximum image found at umL {max_image.nx * r.grid_size_mm.x}, {max_image.ny * r.grid_size_mm.y}"
|
||
)
|
||
com = CenterOfMassModel(n_x=max_image.nx, n_y=max_image.ny, max_image=max_image.number)
|
||
com_mm = com.get_com_mm(r)
|
||
grid_mm_x = com_mm.x
|
||
grid_mm_y = com_mm.y
|
||
|
||
logger.debug(f"Grid coordinates in mm: x={grid_mm_x}, y={grid_mm_y}")
|
||
return com
|
||
else:
|
||
return None
|
||
|
||
|
||
def rebuild_array_from_scan_results(
|
||
scan_results: list,
|
||
value_field: str,
|
||
array_shape: tuple | None = None,
|
||
nx_field: str = "nx",
|
||
ny_field: str = "ny",
|
||
default_value: float = 0.0,
|
||
threshold: float | None = None,
|
||
condition_func: Callable | None = None,
|
||
apply_filter_before: bool = True,
|
||
) -> np.ndarray:
|
||
positions = []
|
||
values = []
|
||
|
||
for result in scan_results:
|
||
nx = getattr(result, nx_field)
|
||
ny = getattr(result, ny_field)
|
||
value = getattr(result, value_field)
|
||
|
||
# Skip if coordinates are None
|
||
if nx is None or ny is None:
|
||
continue
|
||
|
||
positions.append((int(nx), int(ny))) # Note: (row, col) = (ny, nx)
|
||
if not value:
|
||
value = 0.0
|
||
values.append(float(value))
|
||
|
||
if not positions:
|
||
logger.error("No valid positions found in scan results")
|
||
raise ValueError("No valid positions found in scan results")
|
||
|
||
# Determine array shape
|
||
if array_shape is None:
|
||
max_row = max(pos[0] for pos in positions)
|
||
max_col = max(pos[1] for pos in positions)
|
||
array_shape = (max_row + 1, max_col + 1)
|
||
|
||
# Initialize array with default values
|
||
result_array = np.full(array_shape, default_value, dtype=float)
|
||
|
||
# Apply pre-filtering if requested
|
||
if apply_filter_before:
|
||
filtered_data = []
|
||
for pos, val in zip(positions, values):
|
||
keep_value = True
|
||
|
||
# Apply threshold filter
|
||
if threshold is not None and val < threshold:
|
||
keep_value = False
|
||
|
||
# Apply custom condition
|
||
if condition_func is not None and not condition_func(val):
|
||
keep_value = False
|
||
|
||
if keep_value:
|
||
filtered_data.append((pos, val))
|
||
else:
|
||
filtered_data.append((pos, 0.0))
|
||
|
||
# Fill array with filtered values
|
||
for pos, val in filtered_data:
|
||
if 0 <= pos[0] < array_shape[0] and 0 <= pos[1] < array_shape[1]:
|
||
result_array[pos[0], pos[1]] = val
|
||
else:
|
||
# Fill array first, then apply filters
|
||
for pos, val in zip(positions, values):
|
||
if 0 <= pos[0] < array_shape[0] and 0 <= pos[1] < array_shape[1]:
|
||
result_array[pos[0], pos[1]] = val
|
||
|
||
# Apply post-filtering
|
||
if threshold is not None:
|
||
result_array[result_array < threshold] = 0.0
|
||
|
||
if condition_func is not None:
|
||
mask = np.vectorize(condition_func)(result_array)
|
||
result_array[~mask] = 0.0
|
||
|
||
return result_array
|
||
|
||
|
||
def create_quality_filtered_array(
|
||
scan_results: list,
|
||
value_field: str,
|
||
min_spots: int | None = None,
|
||
min_efficiency: float | None = 1.0,
|
||
min_background: float | None = None,
|
||
exclude_ice: bool | None = True,
|
||
min_low_res_spots: float | None = 10.0,
|
||
**kwargs,
|
||
) -> np.ndarray:
|
||
"""
|
||
Create array with comprehensive quality filtering
|
||
"""
|
||
|
||
def quality_condition(
|
||
result,
|
||
max_spots_low_res,
|
||
min_background,
|
||
min_spots,
|
||
min_efficiency,
|
||
min_low_res_spots,
|
||
max_filter: float = 0.3,
|
||
):
|
||
if exclude_ice and (result.spots_ice / max(result.spots_low_res, 1.0)) == 1.0:
|
||
# print(f"all ice for {result.number}")
|
||
return False
|
||
if min_low_res_spots and result.spots_low_res < min_low_res_spots:
|
||
return False
|
||
if result.spots_low_res < (max_spots_low_res * max_filter):
|
||
return False
|
||
if exclude_ice and result.spots_ice > result.spots * 0.8: # More than 50% ice
|
||
# print(f"more than 80% ice for {result.number}")
|
||
return False
|
||
if result.index:
|
||
# print(f"index is True for {result.number}")
|
||
return True
|
||
if min_spots and result.spots < min_spots:
|
||
# print(f"{result.spots} is less than {min_spots} for {result.number}")
|
||
return False
|
||
if result.spots_low_res < min_background:
|
||
return False
|
||
|
||
return not result.efficiency < min_efficiency
|
||
|
||
# Filter results first
|
||
filtered_results = []
|
||
|
||
if min_spots is None:
|
||
min_spots = min(
|
||
(result.spots for result in scan_results if result.spots is not None), default=1
|
||
)
|
||
if min_low_res_spots is None:
|
||
min_low_res_spots = min(
|
||
(result.spots_low_res for result in scan_results if result.spots_low_res is not None),
|
||
default=1,
|
||
)
|
||
if min_background is None:
|
||
min_background = min(
|
||
(result.bkg for result in scan_results if result.bkg is not None), default=1
|
||
)
|
||
if min_efficiency is None:
|
||
min_efficiency = 1.0
|
||
max_spots_low_res = max(
|
||
(result.spots_low_res for result in scan_results if result.spots is not None), default=1
|
||
)
|
||
|
||
for result in scan_results:
|
||
if result.nx is not None and result.ny is not None:
|
||
if quality_condition(
|
||
result,
|
||
max_spots_low_res=max_spots_low_res,
|
||
min_background=min_background,
|
||
min_spots=min_spots,
|
||
min_efficiency=min_efficiency,
|
||
min_low_res_spots=min_low_res_spots,
|
||
):
|
||
filtered_results.append(result)
|
||
else:
|
||
# Create a copy with zero value for filtered positions
|
||
import copy
|
||
|
||
zero_result = copy.copy(result)
|
||
setattr(zero_result, value_field, 0)
|
||
filtered_results.append(zero_result)
|
||
|
||
return rebuild_array_from_scan_results(filtered_results, value_field, **kwargs)
|
||
|
||
|
||
def get_xtal_size(crystal_size, result_array, r: RasterGridRequest):
|
||
# Optional: get bounding box of the largest object
|
||
try:
|
||
labeled_array, num_objects = ndimage.label(result_array)
|
||
areas = ndimage.sum(
|
||
np.ones_like(result_array, dtype=np.int32),
|
||
labeled_array,
|
||
index=range(1, num_objects + 1),
|
||
)
|
||
largest_idx = int(np.argmax(areas)) + 1 # +1 because labels start at 1
|
||
largest_area = int(areas[largest_idx - 1])
|
||
logger.info(f"Largest object label: {largest_idx}, area (px): {largest_area}")
|
||
|
||
object_mask = labeled_array == largest_idx
|
||
rows = np.any(object_mask, axis=1)
|
||
cols = np.any(object_mask, axis=0)
|
||
row_min, row_max = np.where(rows)[0][[0, -1]]
|
||
col_min, col_max = np.where(cols)[0][[0, -1]]
|
||
logger.info(f"Largest bbox: width={col_max - col_min}, height={row_max - row_min}")
|
||
logger.info(
|
||
f"Largest bbox: width={(col_max - col_min) * r.grid_size_mm.x}, y={(row_max - row_min) * r.grid_size_mm.y}"
|
||
)
|
||
if r.n_x == 1:
|
||
logger.debug("calculating z")
|
||
crystal_size = CrystalSize(
|
||
x=crystal_size.x, y=crystal_size.y, z=(col_max - col_min) * r.grid_size_mm.y * 1000
|
||
)
|
||
logger.info(f"Crystal Size: x={crystal_size.x}, y={crystal_size.y}, z={crystal_size.z}")
|
||
|
||
else:
|
||
logger.debug("calculating x and y")
|
||
crystal_size = CrystalSize(
|
||
x=(row_max - row_min) * r.grid_size_mm.x * 1000,
|
||
y=(col_max - col_min) * r.grid_size_mm.y * 1000,
|
||
z=crystal_size.z,
|
||
)
|
||
logger.info(f"Crystal Size: x={crystal_size.x}, y={crystal_size.y}, z={crystal_size.z}")
|
||
except ValueError as e:
|
||
logger.error(f"error calculating xtal size: {e}")
|
||
crystal_size = CrystalSize(x=0, y=0, z=0)
|
||
|
||
return crystal_size
|
||
|
||
|
||
def get_best_b_factor(result_list: list):
|
||
if not result_list:
|
||
return None
|
||
best_b_factor = min(
|
||
(img for img in result_list if img.b is not None), key=lambda img: img.b, default=None
|
||
)
|
||
if best_b_factor is None:
|
||
return None
|
||
logger.info(f"Best b: {best_b_factor.b}")
|
||
return best_b_factor.b
|
||
|
||
|
||
def get_best_res(result_list: list):
|
||
if not result_list:
|
||
return None
|
||
best_res = min(
|
||
(img for img in result_list if img.res is not None), key=lambda img: img.res, default=None
|
||
)
|
||
if best_res is None:
|
||
return None
|
||
logger.info(f"Best res: {best_res.res}")
|
||
return best_res.res
|
||
|
||
|
||
def com_nan_check(com):
|
||
return not (np.isnan(com.n_x) or np.isnan(com.n_y))
|
||
|
||
|
||
def get_result_list_from_com(images, com: CenterOfMassModel):
|
||
if not com_nan_check(com):
|
||
return None
|
||
cx, cy = com.n_x, com.n_y
|
||
start_x, end_x = round(cx - 1), round(cx + 1)
|
||
start_y, end_y = round(cy - 1), round(cy + 1)
|
||
logger.info(f"range x {start_x} {end_x}, y {start_y} {end_y}")
|
||
result_list = [
|
||
img for img in images if start_x <= img.nx <= end_x and start_y <= img.ny <= end_y
|
||
]
|
||
return result_list
|
||
|
||
|
||
def get_com_image_number(com, images):
|
||
for image in images:
|
||
if image.nx == round(com[0]) and image.ny == round(com[1]):
|
||
logger.info(f"com found for image: {image.number}")
|
||
return image.number
|
||
return None
|
||
|
||
|
||
def _to_com_model(coords, images, label: str) -> CenterOfMassModel | None:
|
||
"""Validate (n_x, n_y) grid coords and wrap them in a CenterOfMassModel.
|
||
|
||
Shared by raster_centre_of_mass and raster_highest_score, which differ only
|
||
in how they pick the target cell.
|
||
"""
|
||
if coords and not np.isnan(coords[0]) and not np.isnan(coords[1]):
|
||
logger.info(f"{label}: {coords}")
|
||
max_image = get_com_image_number(coords, images)
|
||
return CenterOfMassModel(n_x=coords[0], n_y=coords[1], max_image=max_image)
|
||
elif coords and np.isnan(coords[0]) and np.isnan(coords[1]):
|
||
logger.warning(f"{label} is nan: {coords[0]}, {coords[1]}")
|
||
return None
|
||
else:
|
||
logger.warning(f"No valid {label} found")
|
||
return None
|
||
|
||
|
||
def raster_centre_of_mass(result_array, images) -> CenterOfMassModel | None:
|
||
# grid_mm_x and grid_mm_y are relative to the top left corner of raster grid
|
||
com = ndimage.center_of_mass(result_array)
|
||
return _to_com_model(com, images, "Center of mass")
|
||
|
||
|
||
def raster_highest_score(images, min_low_res_spots: float = 10.0) -> CenterOfMassModel | None:
|
||
"""Target the grid cell with the highest crystal score.
|
||
|
||
If the whole grid is too weak to hold a crystal — max spots_low_res below
|
||
min_low_res_spots — target the geometric centre of the grid instead of
|
||
collecting at a noisy cell, so a 'nothing here' result is centred and
|
||
deliberate rather than random noise.
|
||
|
||
spots_indexed is deliberately not part of the guard: indexing was dropped
|
||
from the crystal score (w_indexed=0.00), so a crystal that diffracts but
|
||
fails to index must still be targeted, not routed to centre.
|
||
"""
|
||
score_array = compute_crystal_score_array(images)
|
||
max_low_res = max((getattr(img, "spots_low_res", 0) or 0 for img in images), default=0)
|
||
if max_low_res < min_low_res_spots:
|
||
n_nx, n_ny = score_array.shape
|
||
logger.info(
|
||
f"No crystal in loop (max spots_low_res={max_low_res} < {min_low_res_spots}); "
|
||
f"targeting grid centre of {n_nx}x{n_ny} grid"
|
||
)
|
||
# get_com_mm maps index i -> (i+0.5)*step, so index (N-1)/2 is the true
|
||
# geometric centre of the grid for both odd and even N (and N==1).
|
||
return CenterOfMassModel(n_x=(n_nx - 1) / 2.0, n_y=(n_ny - 1) / 2.0)
|
||
return _to_com_model(_max_cell(score_array), images, "Highest score")
|
||
|
||
|
||
def has_sufficient_low_res_spots(result_array: np.ndarray, min_spots_low_res: float) -> bool:
|
||
if result_array is None or result_array.size == 0:
|
||
logger.warning("Empty result_array passed to has_sufficient_low_res_spots")
|
||
return False
|
||
|
||
max_val = float(np.nanmax(result_array))
|
||
logger.info(f"Max spots_low_res in raster result_array: {max_val}")
|
||
|
||
if max_val < min_spots_low_res:
|
||
logger.info(f"Raster rejected: max spots_low_res {max_val} < required {min_spots_low_res}")
|
||
return False
|
||
|
||
return True
|
||
|
||
|
||
def compute_crystal_score_array(
|
||
scan_results: list, w_bkg: float = 0.25, w_low_res: float = 0.75, w_indexed: float = 0.00
|
||
) -> np.ndarray:
|
||
"""Combine bkg (25%), spots_low_res (75%), and spots_indexed (00%) into a 0–100 score.
|
||
|
||
Each field is min-max normalised to [0, 100] within the grid before weighting,
|
||
so the final score is the probability (0–100) that a pixel belongs to a crystal.
|
||
"""
|
||
arr_bkg = rebuild_array_from_scan_results(scan_results, "bkg")
|
||
arr_low = rebuild_array_from_scan_results(scan_results, "spots_low_res")
|
||
arr_idx = rebuild_array_from_scan_results(scan_results, "spots_indexed")
|
||
|
||
def _norm(a: np.ndarray) -> np.ndarray:
|
||
mn, mx = float(a.min()), float(a.max())
|
||
if mx == mn:
|
||
return np.zeros_like(a, dtype=float)
|
||
return (a - mn) / (mx - mn) * 100.0
|
||
|
||
return w_bkg * _norm(arr_bkg) + w_low_res * _norm(arr_low) + w_indexed * _norm(arr_idx)
|
||
|
||
|
||
def _max_cell(arr: np.ndarray) -> tuple[int, int]:
|
||
"""Return (nx, ny) of the cell with the highest value."""
|
||
idx = np.unravel_index(np.argmax(arr), arr.shape)
|
||
return int(idx[0]), int(idx[1])
|
||
|
||
|
||
def _draw_panel(
|
||
ax,
|
||
arr: np.ndarray,
|
||
mask: np.ndarray,
|
||
label: str,
|
||
threshold: float,
|
||
grid_size_mm: tuple[float, float] | None = None,
|
||
cbar_label: str = "spots_low_res",
|
||
) -> None:
|
||
"""Shared helper: heatmap + crystal contour + max-cell square on one Axes.
|
||
|
||
If grid_size_mm=(step_x_mm, step_y_mm) is provided, the crystal size in µm
|
||
is computed via get_xtal_size and shown in the panel title.
|
||
"""
|
||
import matplotlib.patches as mpatches
|
||
import matplotlib.pyplot as plt
|
||
|
||
from aarecommon.coordinate import Coordinate
|
||
|
||
n_nx, n_ny = arr.shape
|
||
max_nx, max_ny = _max_cell(arr)
|
||
n_cells = int(mask.sum())
|
||
|
||
im = ax.imshow(
|
||
arr.T,
|
||
origin="lower",
|
||
cmap="viridis",
|
||
aspect="equal",
|
||
extent=[-0.5, n_nx - 0.5, -0.5, n_ny - 0.5],
|
||
)
|
||
plt.colorbar(im, ax=ax, label=cbar_label, fraction=0.046, pad=0.04)
|
||
ax.contour(
|
||
np.arange(n_nx),
|
||
np.arange(n_ny),
|
||
mask.T.astype(float),
|
||
levels=[0.5],
|
||
colors="cyan",
|
||
linewidths=1.8,
|
||
)
|
||
rect = mpatches.Rectangle(
|
||
(max_nx - 0.5, max_ny - 0.5), 1, 1, linewidth=2, edgecolor="red", facecolor="none"
|
||
)
|
||
ax.add_patch(rect)
|
||
|
||
size_line = ""
|
||
if grid_size_mm is not None and n_cells > 0:
|
||
r = RasterGridRequest(
|
||
exp_time_s=0.0,
|
||
n_x=n_nx,
|
||
n_y=n_ny,
|
||
grid_size_mm=Coordinate(x=grid_size_mm[0], y=grid_size_mm[1]),
|
||
smargon_top_left=None,
|
||
)
|
||
xtal_size = get_xtal_size(CrystalSize(x=0, y=0, z=0), mask.astype(float), r)
|
||
size_line = f"\n{xtal_size.x:.0f}×{xtal_size.y:.0f} µm"
|
||
|
||
ax.set_title(f"{label}\nthresh≈{threshold:.1f} | {n_cells} cells{size_line}", fontsize=8)
|
||
ax.set_xlabel("nx", fontsize=7)
|
||
ax.set_ylabel("ny", fontsize=7)
|
||
ax.tick_params(labelsize=6)
|
||
|
||
|
||
# ── Clustering / thresholding methods ─────────────────────────────────────────
|
||
|
||
|
||
def crystal_mask_corner_background(arr: np.ndarray) -> tuple[np.ndarray, float]:
|
||
"""Threshold = max(far-corner value, floor=10).
|
||
|
||
Simple, parameter-free. Returns too many cells when corners are zero
|
||
because any nonzero value passes.
|
||
"""
|
||
n_nx, n_ny = arr.shape
|
||
corners = [arr[0, 0], arr[n_nx - 1, 0], arr[0, n_ny - 1], arr[n_nx - 1, n_ny - 1]]
|
||
threshold = max(float(np.max(corners)), 10.0)
|
||
return arr > threshold, threshold
|
||
|
||
|
||
def crystal_mask_otsu_nonzero(arr: np.ndarray) -> tuple[np.ndarray, float]:
|
||
"""Otsu threshold computed only on the nonzero values.
|
||
|
||
Finds the natural gap in the signal distribution.
|
||
Ignores the large mass of background zeros so the threshold is
|
||
placed within the diffraction-signal population.
|
||
"""
|
||
from skimage.filters import threshold_otsu
|
||
|
||
nonzero = arr[arr > 0]
|
||
if nonzero.size == 0:
|
||
return np.zeros_like(arr, dtype=bool), 0.0
|
||
threshold = float(threshold_otsu(nonzero))
|
||
return arr > threshold, threshold
|
||
|
||
|
||
def crystal_mask_mean_sigma(arr: np.ndarray, n_sigma: float = 0.5) -> tuple[np.ndarray, float]:
|
||
"""Threshold = mean + n_sigma * std of nonzero values.
|
||
|
||
n_sigma=0.5 keeps cells within ~1/2 std above average signal.
|
||
Raise n_sigma to tighten the region around the strongest-diffracting core.
|
||
"""
|
||
nonzero = arr[arr > 0]
|
||
if nonzero.size == 0:
|
||
return np.zeros_like(arr, dtype=bool), 0.0
|
||
threshold = float(nonzero.mean() + n_sigma * nonzero.std())
|
||
return arr > threshold, threshold
|
||
|
||
|
||
def crystal_mask_signal_percentile(
|
||
arr: np.ndarray, percentile: float = 60.0
|
||
) -> tuple[np.ndarray, float]:
|
||
"""Threshold = given percentile of nonzero values.
|
||
|
||
percentile=60 keeps the top 40 % of signal cells; raise to sharpen.
|
||
Unlike mean+sigma this is robust to heavy-tailed distributions.
|
||
"""
|
||
nonzero = arr[arr > 0]
|
||
if nonzero.size == 0:
|
||
return np.zeros_like(arr, dtype=bool), 0.0
|
||
threshold = float(np.percentile(nonzero, percentile))
|
||
return arr > threshold, threshold
|
||
|
||
|
||
def crystal_mask_dbscan(
|
||
arr: np.ndarray, min_signal: float = 10.0, eps: float = 1.5, min_samples: int = 3
|
||
) -> tuple[np.ndarray, float]:
|
||
"""DBSCAN spatial clustering on cells with signal > min_signal.
|
||
|
||
Groups adjacent diffracting cells into clusters; the largest cluster
|
||
(by total signal weight) is labelled as the crystal.
|
||
eps=1.5 connects cells that are one grid step apart (including diagonal).
|
||
"""
|
||
from sklearn.cluster import DBSCAN
|
||
|
||
ys, xs = np.where(arr > min_signal)
|
||
if len(ys) == 0:
|
||
return np.zeros_like(arr, dtype=bool), min_signal
|
||
|
||
coords = np.column_stack([ys, xs]).astype(float)
|
||
labels = DBSCAN(eps=eps, min_samples=min_samples).fit_predict(coords)
|
||
|
||
best_label, best_weight = -1, -1.0
|
||
for lbl in set(labels):
|
||
if lbl == -1:
|
||
continue
|
||
weight = float(arr[ys[labels == lbl], xs[labels == lbl]].sum())
|
||
if weight > best_weight:
|
||
best_label, best_weight = lbl, weight
|
||
|
||
mask = np.zeros_like(arr, dtype=bool)
|
||
if best_label != -1:
|
||
sel = labels == best_label
|
||
mask[ys[sel], xs[sel]] = True
|
||
return mask, min_signal
|
||
|
||
|
||
def crystal_mask_kmeans(arr: np.ndarray, n_clusters: int = 3) -> tuple[np.ndarray, float]:
|
||
"""K-means on signal values (1-D feature).
|
||
|
||
Splits cells into n_clusters groups; the cluster with the highest
|
||
centroid is the crystal. n_clusters=3 separates background / fringe / crystal.
|
||
"""
|
||
from sklearn.cluster import KMeans
|
||
|
||
flat = arr.flatten().reshape(-1, 1)
|
||
km = KMeans(n_clusters=n_clusters, random_state=0, n_init="auto").fit(flat)
|
||
crystal_label = int(np.argmax(km.cluster_centers_.flatten()))
|
||
threshold = float(sorted(km.cluster_centers_.flatten())[-2])
|
||
mask = km.labels_.reshape(arr.shape) == crystal_label
|
||
return mask, threshold
|
||
|
||
|
||
CRYSTAL_METHOD_MAP: dict = {
|
||
"corner": ("Corner background", crystal_mask_corner_background),
|
||
"otsu": ("Otsu (nonzero)", crystal_mask_otsu_nonzero),
|
||
"mean_sigma": ("Mean + 0.5σ", lambda arr: crystal_mask_mean_sigma(arr, n_sigma=0.5)),
|
||
"percentile": (
|
||
"Top-40% signal",
|
||
lambda arr: crystal_mask_signal_percentile(arr, percentile=60.0),
|
||
),
|
||
"dbscan": ("DBSCAN (spatial)", crystal_mask_dbscan),
|
||
"kmeans": ("K-means (k=3)", crystal_mask_kmeans),
|
||
}
|
||
|
||
|
||
def compare_crystal_methods(
|
||
results: list, arr: np.ndarray | None = None, grid_size_mm: tuple[float, float] | None = None
|
||
) -> None:
|
||
"""Plot each crystal-detection method side-by-side for visual comparison.
|
||
|
||
Cyan contour = detected crystal region.
|
||
Red square = cell with maximum spots_low_res.
|
||
grid_size_mm = (step_x_mm, step_y_mm) — when provided, crystal size in µm
|
||
is calculated via get_xtal_size and shown in each panel title.
|
||
"""
|
||
import matplotlib.pyplot as plt
|
||
|
||
if arr is None:
|
||
arr = rebuild_array_from_scan_results(results, "spots_low_res")
|
||
|
||
methods = [
|
||
("Corner background", *crystal_mask_corner_background(arr)),
|
||
("Otsu (nonzero)", *crystal_mask_otsu_nonzero(arr)),
|
||
("Mean + 0.5σ", *crystal_mask_mean_sigma(arr, n_sigma=0.5)),
|
||
("Top-40% signal", *crystal_mask_signal_percentile(arr, percentile=60.0)),
|
||
("DBSCAN (spatial)", *crystal_mask_dbscan(arr)),
|
||
("K-means (k=3)", *crystal_mask_kmeans(arr)),
|
||
]
|
||
|
||
fig, axes = plt.subplots(2, 3, figsize=(14, 9))
|
||
for ax, (label, mask, threshold) in zip(axes.flat, methods):
|
||
_draw_panel(ax, arr, mask, label, threshold, grid_size_mm=grid_size_mm)
|
||
|
||
fig.suptitle(
|
||
"Crystal region detection — method comparison\n"
|
||
"cyan = crystal outline | red square = max spots_low_res cell",
|
||
fontsize=10,
|
||
)
|
||
plt.tight_layout()
|
||
plt.show()
|
||
|
||
|
||
def compare_crystal_methods_scored(
|
||
results: list,
|
||
score_arr: np.ndarray | None = None,
|
||
method: str | None = None,
|
||
grid_size_mm: tuple[float, float] | None = None,
|
||
w_bkg: float = 0.20,
|
||
w_low_res: float = 0.60,
|
||
w_indexed: float = 0.20,
|
||
) -> None:
|
||
"""Plot each crystal-detection method applied to the composite 0–100 crystal score.
|
||
|
||
The score combines bkg (w_bkg), spots_low_res (w_low_res), and spots_indexed
|
||
(w_indexed), each min-max normalised. All six methods are shown side-by-side.
|
||
|
||
Args:
|
||
results: scan result list (used to build score if score_arr is None)
|
||
score_arr: pre-computed score array; computed from results when None
|
||
method: if given (one of 'corner','otsu','mean_sigma','percentile',
|
||
'dbscan','kmeans'), that panel is highlighted with a yellow border
|
||
grid_size_mm: (step_x_mm, step_y_mm) for crystal-size annotation
|
||
w_bkg, w_low_res, w_indexed: composite score weights (must sum to 1.0)
|
||
"""
|
||
import matplotlib.pyplot as plt
|
||
|
||
if score_arr is None:
|
||
score_arr = compute_crystal_score_array(
|
||
results, w_bkg=w_bkg, w_low_res=w_low_res, w_indexed=w_indexed
|
||
)
|
||
|
||
method_entries = [(name, label, fn) for name, (label, fn) in CRYSTAL_METHOD_MAP.items()]
|
||
|
||
fig, axes = plt.subplots(2, 3, figsize=(15, 9))
|
||
for ax, (name, label, fn) in zip(axes.flat, method_entries):
|
||
mask, threshold = fn(score_arr)
|
||
_draw_panel(
|
||
ax,
|
||
score_arr,
|
||
mask,
|
||
label,
|
||
threshold,
|
||
grid_size_mm=grid_size_mm,
|
||
cbar_label="crystal score (0–100)",
|
||
)
|
||
if method and name == method:
|
||
for spine in ax.spines.values():
|
||
spine.set_edgecolor("yellow")
|
||
spine.set_linewidth(3)
|
||
|
||
weight_note = (
|
||
f"bkg×{w_bkg:.0%} + spots_low_res×{w_low_res:.0%} + spots_indexed×{w_indexed:.0%}"
|
||
)
|
||
fig.suptitle(
|
||
f"Crystal region detection on composite score [{weight_note}]\n"
|
||
"cyan = crystal outline | red square = max-score cell | score range 0–100",
|
||
fontsize=10,
|
||
)
|
||
plt.tight_layout()
|
||
plt.show()
|
||
|
||
|
||
if __name__ == "__main__":
|
||
from types import SimpleNamespace as _I
|
||
|
||
def _cell(nx, ny, low, idx, bkg=0):
|
||
return _I(nx=nx, ny=ny, number=nx * 3 + ny, spots_low_res=low, spots_indexed=idx, bkg=bkg)
|
||
|
||
# Real crystal at (1,1) -> best cell returned.
|
||
strong = [_cell(x, y, 0, 0) for x in range(3) for y in range(3)]
|
||
strong[4] = _cell(1, 1, 50, 5, bkg=100)
|
||
com = raster_highest_score(strong)
|
||
assert (round(com.n_x), round(com.n_y)) == (1, 1), com
|
||
|
||
# All noise (low-res < 10, nothing indexed) -> centre of 3x3 grid == (1,1).
|
||
noise = [_cell(x, y, 3, 0) for x in range(3) for y in range(3)]
|
||
com = raster_highest_score(noise)
|
||
assert (com.n_x, com.n_y) == (1.0, 1.0), com
|
||
|
||
# Strong low-res but nothing indexed anywhere -> still targets the best
|
||
# cell: indexing is dropped from the score (w_indexed=0.00) and must not
|
||
# veto a diffracting crystal.
|
||
no_index = [_cell(x, y, 12, 0) for x in range(3) for y in range(3)]
|
||
no_index[6] = _cell(2, 0, 40, 0)
|
||
com = raster_highest_score(no_index)
|
||
assert (round(com.n_x), round(com.n_y)) == (2, 0), com
|
||
|
||
print("find_xtal self-check passed")
|