Merge branch 'daq_raster_predicition'

# Conflicts:
#	daq/src/aaredaq/daq.py
#	gui/src/aaregui/panels/rotation_data_collection.py
This commit is contained in:
2025-10-20 12:48:47 +02:00
11 changed files with 468 additions and 404 deletions
Binary file not shown.
+43 -314
View File
@@ -6,6 +6,7 @@ from datetime import datetime
from math import ceil
from typing import List, Tuple, Optional, Callable, Dict
import secrets
import os
import cv2
import numpy as np
@@ -23,12 +24,14 @@ from aaredaqlib.beamline import MXBeamline
from aaredaqlib.coordinate import Coordinate, SmargonCoordinate
from aaredaqlib.diffraction_geometry import DiffractionGeometry
from aaredaqlib.logger_config import setup_logger
from aaredaqlib.find_xtal import raster_centre_of_mass, create_quality_filtered_array, identify_crystal_raster, \
get_result_list_from_com, get_best_b_factor, get_best_res, get_xtal_size
from aaredaqlib.models import (
SampleShortInfo,
PuckLoadedInfo,
SampleShortInfoList,
DAQStatusModel, BeamlineStatus, SessionStatus, SampleCameraSettings, AutofocusSettings, BoundingBoxModel,
ZoomModeEnum, CrystalSize)
ZoomModeEnum, CrystalSize, SimpleScanParameters)
from aaredaqlib.raster_grid import RasterGridRequest, CompletedRasterGrid, CompletedRasterGridElem
from aaredaqlib.rotation_scan import RotationScanRequest, CompletedRotationScan
from aaredaqlib.sample_geometry import SampleGeometryModel
@@ -64,6 +67,7 @@ class AareDAQ:
self.crystal_size = CrystalSize(x=0,y=0,z=0)
self.last_best_b_factor = None
self.last_best_res = None
self.auto_params = SimpleScanParameters()
@property
def state(self) -> BeamlineStateEnum:
@@ -88,6 +92,14 @@ class AareDAQ:
self.last_time = end - start
@property
def smart_params(self) -> SimpleScanParameters:
return self.auto_params
@smart_params.setter
def smart_params(self, params: SimpleScanParameters):
self.auto_params = params
@property
def omega(self) -> float:
return self.__devs.aerotech.omega
@@ -369,202 +381,22 @@ class AareDAQ:
def list_loaded_pucks(self) -> List[PuckLoadedInfo]:
return self.__devs.tell.get_detected_pucks()
def identify_crystal_raster(self, result, r: RasterGridRequest):
images = result.images
if images and any(img.spots > 0 for img in images):
indexed_images = [img for img in images if img.index and img.spots_low_res > 4 and img.bkg > 4.5]
if indexed_images:
#indexed_images = [img for img in indexed_images if img.spots_indexed > 10]
images = indexed_images
filtered_images = [img for img in images
if img.spots_ice is not None and img.spots_low_res > 4 and (
img.spots_ice / img.spots_low_res) < 5.0
and (img.spots_ice / img.spots_low_res) != 1]
if indexed_images:
print(f"Find image by maximum number of spots indexed")
max_image = max(images, key=lambda img: img.spots_indexed)
max_spots = max_image.spots_indexed
max_images = [img for img in images if img.spots_indexed == max_spots]
max_image = max_images[len(max_images) // 2]
print(f"Image with maximum spots_low_res: {max_image}")
print(f"Maximum spots_indexed value: {max_image.spots_indexed}")
print(f"Maximum spots_low_res value: {max_image.spots_low_res}")
else:
print(f"Find image by maximum number of low resolution spots")
max_image = max(images, key=lambda img: img.spots_low_res)
print(f"Image with maximum spots_low_res: {max_image}")
print(f"Maximum spots_low_res value: {max_image.spots_low_res}")
grid_mm_x = max_image.nx * r.grid_size_mm.x
grid_mm_y = max_image.ny * r.grid_size_mm.y
print(f"Grid coordinates in mm: ({grid_mm_x}, {grid_mm_y})")
delta_mm = self.sample_geometry.smargon_nudge(Coordinate(x=grid_mm_x, y=grid_mm_y))
return delta_mm
else:
return None
def rebuild_array_from_scan_results(self,
scan_results: List,
value_field: str,
array_shape: Optional[tuple] = None,
nx_field: str = 'nx',
ny_field: str = 'ny',
default_value: float = 0.0,
threshold: Optional[float] = None,
condition_func: Optional[Callable] = None,
apply_filter_before: bool = True
) -> np.ndarray:
# Extract coordinates and values
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:
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(self,
scan_results: List,
value_field: str,
min_spots: Optional[int] = None,
min_efficiency: Optional[float] = 1.0,
min_background: Optional[float] = None,
exclude_ice: Optional[bool] = True,
**kwargs
) -> np.ndarray:
"""
Create array with comprehensive quality filtering
"""
def quality_condition(result, min_bkg, min_spots, min_efficiency):
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 (result.spots_ice / result.spots_low_res) > 5.0:
# 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 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
if result.efficiency < min_efficiency:
return False
return True
# 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_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
for result in scan_results:
if result.nx is not None and result.ny is not None:
if quality_condition(result, min_bkg=min_background, min_spots=min_spots,
min_efficiency=min_efficiency):
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 self.rebuild_array_from_scan_results(filtered_results, value_field, **kwargs)
def __auto_center(self, grid: RasterGridRequest) -> CompletedRasterGrid | None:
sample = self.sample
if sample is None:
raise Exception("Sample must be mounted to auto center")
old_prefix = grid.file_prefix
geom = self.sample_geometry
r = self.__ml_bounding_box(sample.db_id, f"ml_{geom.omega_deg:.2f}deg")
if r is None:
self.__devs.aerotech.move(geom.omega_deg + 90.0, wait=True)
time.sleep(0.2)
r = self.__ml_bounding_box(sample.db_id, f"ml_{geom.omega_deg + 90.0:.2f}deg")
if r is not None:
geom = self.sample_geometry
grid.smargon = r.smargon
@@ -590,132 +422,6 @@ class AareDAQ:
else:
return None
def raster_centre_of_mass(self, images, r:RasterGridRequest, result):
# if any(img.index for img in images):
# print("COM by indexed spots")
# result_array= self.create_quality_filtered_array(images, 'spots_indexed', min_spots=None,
# min_efficiency=1.0, min_background=None)
# else:
print("COM by low res spots")
result_array = self.create_quality_filtered_array(images, 'spots_low_res', min_spots=None,
min_efficiency=1.0, min_background=None)
print('horizontal scan')
com = ndimage.center_of_mass(result_array)
print(f"Center of mass: {com}")
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])
print(f"Largest object label: {largest_idx}, area (px): {largest_area}")
# Optional: get bounding box of largest object
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]]
print(f"Largest bbox: width={col_max - col_min}, height={row_max - row_min}")
print(
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:
if self.crystal_size is None:
self.crystal_size = CrystalSize(x=0,y=0,z=0)
crystal_size = self.crystal_size
crystal_size = CrystalSize(x=crystal_size.x, y=crystal_size.y,
z=(col_max - col_min) * r.grid_size_mm.y * 1000)
else:
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=0)
except Exception as e:
crystal_size = CrystalSize(x=0, y=0, z=0)
self.crystal_size = crystal_size
if r.n_x == 1 and (np.isnan(com[1]) or np.isnan(com[0])):
print('vertical scan')
try:
max_image = max(images, key=lambda img: img.spots_low_res)
com = (0, max_image.ny)
except:
print("no spots")
if np.isnan(com[1]) or np.isnan(com[0]):
print("Center of mass is nan")
com = None
grid_mm_x = None
grid_mm_y = None
best_res = None
best_b_factor = None
else:
grid_mm_x = com[0] * r.grid_size_mm.x
grid_mm_y = com[1] * r.grid_size_mm.y
if r.n_x == 1:
grid_mm_x += (0.5 * r.grid_size_mm.x)
cx, cy = com[::-1] # com=(y, x) -> (x, y), rounded once
start_x, end_x = round(cx - 1), round(cx + 1)
start_y, end_y = round(cy - 1), round(cy + 1)
print(f"range x {start_x} {end_x}, y {start_y} {end_y}")
# Collect images in the 3x3 neighborhood around the center
res_list = [img for img in images
if start_x <= img.nx <= end_x and start_y <= img.ny <= end_y]
print(res_list)
# Best by res, skipping None
best_res = min((img for img in res_list if img.res is not None),
key=lambda img: img.res,
default=None)
print(f"Best res: {best_res}")
best_b_factor = min((img for img in res_list if img.b is not None),
key=lambda img: img.b,
default=None)
print(f"Best b: {best_b_factor}")
for image in images:
if image.nx == round(com[0]) and image.ny == round(com[1]):
print(f"com found for image: {image.number}")
try:
if best_res is not None and best_res.res is not None:
# store the numeric resolution on the config/session so it appears in status
print(f"best res: {best_res.res}")
self.last_best_res = float(best_res.res)
else:
print(f'res is None')
self.last_best_res = None
except Exception as e:
print(f'error with last_best_res {e}')
self.last_best_res = None
try:
if best_b_factor is not None and best_b_factor.b is not None:
print(f"best res: {best_res.b}")
self.last_best_b_factor = float(best_b_factor.b)
else:
print(f'last_best_b_factor is None')
self.last_best_b_factor = None
except Exception as e:
print(f'error with last_best_b_factor {e}')
self.last_best_b_factor = None
if com is not None and grid_mm_x is not None and grid_mm_y is not None:
new_delta_mm = self.sample_geometry.smargon_nudge(Coordinate(x=grid_mm_x, y=grid_mm_y))
print(f"new delta mm: {new_delta_mm}, new grid x: {grid_mm_x}, new grid y: {grid_mm_y}")
else:
print(f"using old method as COM is none or nan")
new_delta_mm = self.identify_crystal_raster(result, r)
return new_delta_mm, best_res
def __raster(self, r: RasterGridRequest) -> CompletedRasterGridElem:
max_time = r.exp_time_s * r.n_y * r.n_x + 60
@@ -774,7 +480,28 @@ class AareDAQ:
with open(filename, 'w') as f:
json.dump(output_data, f, indent=2)
print('before centre_of_mass')
new_delta_mm, best_res = self.raster_centre_of_mass(images, r, result)
print("COM by low res spots")
result_array = create_quality_filtered_array(images, 'spots_low_res', min_spots=None,
min_efficiency=1.0, min_background=None, min_low_res_spots=10.0)
self.crystal_size = get_xtal_size(self.crystal_size, result_array, r)
grid_mm_x, grid_mm_y, com = raster_centre_of_mass(result_array, r)
if grid_mm_x is None or grid_mm_y is None:
print(f"using old method as COM is none or nan")
grid_mm_x, grid_mm_y = identify_crystal_raster(result, r)
if r.n_x == 1:
new_delta_mm = self.sample_geometry.smargon_nudge(Coordinate(x=0, y=grid_mm_y))
else:
result_list = get_result_list_from_com(images, com)
print(result_list)
self.last_best_b_factor = get_best_b_factor(result_list)
self.last_best_res = get_best_res(result_list)
print("b_factor: ", self.last_best_b_factor, " best_res: ", self.last_best_res)
new_delta_mm = self.sample_geometry.smargon_nudge(Coordinate(x=grid_mm_x, y=grid_mm_y))
print('after centre_of_mass')
if new_delta_mm is not None:
print(f'{time.ctime()}, moving SMARGON to target new delta mm {r.smargon.sh_mm + new_delta_mm} mm')
@@ -1490,7 +1217,6 @@ class AareDAQ:
def measure(self, sample: SampleShortInfo) -> float:
start = time.perf_counter()
formatted_date = datetime.now().strftime('%Y%m%d')
sample_prefix = "{}/{}/{:02d}/{}".format(
formatted_date,
@@ -1533,6 +1259,8 @@ class AareDAQ:
grid_size_mm=Coordinate(x=geom.beam_size_mm.x * 0.5, y=geom.beam_size_mm.y * 0.5),
omega_deg=0
)):
params = self.smart_params
print(params)
self.__cfg.zoom_mode = ZoomModeEnum.User
self.__devs.samcam_settings = self.__cfg.zoom_settings.get_camera_settings(self.zoom)
self.__rotation( RotationScanRequest(start_omega_deg=0,
@@ -1542,13 +1270,14 @@ class AareDAQ:
incr_omega_deg=0.2,
steps=900,))
else:
print("auto center failed")
self.__aare.axc_failed(sample)
self.zoom = 1
self.__cfg.zoom_mode = ZoomModeEnum.User
self.__devs.samcam_settings = self.__cfg.zoom_settings.get_camera_settings(self.zoom)
self.__cfg.state_busy = False
except Exception:
logger.exception("Exception in measure")
except Exception as e:
print(f"Error in measure: {e}")
self.__aare.sample_failed(sample)
self.__cfg.state_busy = False
end = time.perf_counter()
+7 -1
View File
@@ -8,7 +8,7 @@ import uvicorn
from aaredaqlib.coordinate import SmargonCoordinate, Coordinate
from aaredaqlib.models import SampleShortInfo, DAQStatusModel, BeamlineStateEnum, BeamlineSettingsModel, \
SampleShortInfoList, SessionStatus, SampleCameraSettings, AutofocusSettings, TokenData, \
CryojetSettingsModel
CryojetSettingsModel, SimpleScanParameters
from aaredaqlib.raster_grid import RasterGridRequest, CompletedRasterGrid
from aaredaqlib.rotation_scan import RotationScanRequest, CompletedRotationScan
from aaredaqlib.sample_geometry import SampleGeometryModel
@@ -345,6 +345,12 @@ async def auto(s: SampleShortInfo, token: str = Depends(oauth2_scheme)):
runtime = daq.measure(s)
return f"{runtime:0.3f}"
@app.post("/scan/smart_params")
async def set_smart_params(p: SimpleScanParameters, token: str = Depends(oauth2_scheme)) -> str:
auth.check_jwt_rw(cfg, auth.parse_token(token))
daq.smart_params = p
return "OK"
@app.post("/scan/cancel")
async def cancel(token: str = Depends(oauth2_scheme)):
auth.check_jwt_rw(cfg, auth.parse_token(token))
+2 -3
View File
@@ -143,6 +143,7 @@ def dc2rse(devs: BeamlineDevices, cfg: BeamlineConfig):
def se2sa(devs: BeamlineDevices, cfg: BeamlineConfig):
hub = cfg.settings
cfg.zoom_mode = ZoomModeEnum.User
devs.detector_cover.put(1)
devs.samcam_settings = cfg.zoom_settings.get_camera_settings(devs.zoom)
hub_cryo = cfg.cryojet_settings
_cryo = hub_cryo.cryojet_in_use
@@ -152,7 +153,7 @@ def se2sa(devs: BeamlineDevices, cfg: BeamlineConfig):
devs.aerotech.unlock()
devs.aerotech.set_direct_mode()
print(time.ctime(), " moving reflector to up position")
devs.reflector_up = True
#devs.reflector_up = True
print(time.ctime(), " moving beamstop to up position")
devs.beamstop_stage_up = True
print(time.ctime(), " setting lamp to 2.5")
@@ -167,8 +168,6 @@ def se2sa(devs: BeamlineDevices, cfg: BeamlineConfig):
devs.dtz.move(cfg.dtz, wait=False)
def rse2sa(devs: BeamlineDevices, cfg: BeamlineConfig):
devs.detector_cover.put(2) # Open detector cover
print(time.ctime(), " moving DETECTOR COVER to OPEN position")
if devs.tell.is_in_park():
print(time.ctime(), " moving TELL to COLD")
devs.tell.move_cold(wait=True)
+5 -2
View File
@@ -126,8 +126,11 @@ class JFJochWrapper:
ring_current_mA=s.bl.ring_current_mA,
sample_temperature_K=s.bl.cryojet_K,
total_flux=s.bl.flux_ph_s,
space_group_number=1
#TODO space_group_number and unit_cell???
space_group_number=1,#TODO space_group_number and unit_cell???
poni_rot1_rad = s.diffraction.poni_rot1_rad,
poni_rot2_rad = s.diffraction.poni_rot2_rad,
max_spot_count = 1000,
detect_ice_rings = True
)
self.__api.start_post(dataset_settings=dataset_settings)