May 2025
This commit is contained in:
root
@@ -79,7 +79,7 @@ def edge(filter_name, backgrounds, signals, peakback):
|
||||
sig_deriv -= peakback
|
||||
sig_deriv *= signal.tukey(2048) # just added Nov 20 2024
|
||||
#peak_pos = 1024 - (np.argmax(sig_deriv[500:1500], axis=-1) + 500)
|
||||
peak_pos = 800 - (np.argmax(sig_deriv[400:1600], axis=-1) + 400) # I am grossed out
|
||||
peak_pos = 1024 - (np.argmax(sig_deriv[400:1600], axis=-1) + 400) # I am grossed out
|
||||
peak_amp = np.amax(sig_deriv[400:1600], axis=-1)
|
||||
|
||||
return peak_pos, peak_amp, sig_deriv, sig_uninterp
|
||||
|
||||
@@ -25,6 +25,12 @@ base_pv_names = []
|
||||
all_pv_names = []
|
||||
global_ravg_length = 100
|
||||
ravg_buffers = {}
|
||||
# Configuration for rolling-average of statistics
|
||||
|
||||
# Configuration for N-shot average spectrum
|
||||
global_avg_length = 100
|
||||
avg_buffer = None
|
||||
avg_pv_names = []
|
||||
|
||||
@numba.njit(parallel=False)
|
||||
def get_spectrum(image, background):
|
||||
@@ -55,51 +61,66 @@ def update_PVs(buffer, *pv_names):
|
||||
|
||||
|
||||
def initialize(params):
|
||||
"""Initialize PV names, running-average settings, and launch update thread."""
|
||||
"""Initialize PV names, running-average settings, N-shot average, and launch update thread."""
|
||||
global channel_pv_names, base_pv_names, all_pv_names, global_ravg_length
|
||||
global global_avg_length, avg_buffer, avg_pv_names
|
||||
|
||||
camera = params["camera_name"]
|
||||
e_int = params["e_int_name"]
|
||||
e_axis = params["e_axis_name"]
|
||||
|
||||
# Fit/result PV names
|
||||
center_pv = f"{camera}:FIT-COM"
|
||||
fwhm_pv = f"{camera}:FIT-FWHM"
|
||||
fit_rms_pv = f"{camera}:FIT-RMS"
|
||||
fit_res_pv = f"{camera}:FIT-RES"
|
||||
center_pv = f"{camera}:FIT-COM"
|
||||
fwhm_pv = f"{camera}:FIT-FWHM"
|
||||
fit_rms_pv = f"{camera}:FIT-RMS"
|
||||
fit_res_pv = f"{camera}:FIT-RES"
|
||||
fit_spec_pv = f"{camera}:FIT-SPECTRUM_Y"
|
||||
|
||||
# ROI PVs for dynamic read
|
||||
ymin_pv = f"{camera}:SPC_ROI_YMIN"
|
||||
ymax_pv = f"{camera}:SPC_ROI_YMAX"
|
||||
axis_pv = e_axis
|
||||
ymin_pv = f"{camera}:SPC_ROI_YMIN"
|
||||
ymax_pv = f"{camera}:SPC_ROI_YMAX"
|
||||
axis_pv = e_axis
|
||||
channel_pv_names = [ymin_pv, ymax_pv, axis_pv]
|
||||
|
||||
# Spectrum statistical PV names
|
||||
com_pv = f"{camera}:SPECT-COM"
|
||||
std_pv = f"{camera}:SPECT-RMS"
|
||||
skew_pv = f"{camera}:SPECT-SKEW"
|
||||
iqr_pv = f"{camera}:SPECT-IQR"
|
||||
res_pv = f"{camera}:SPECT-RES" # will use IQR-based calc
|
||||
com_pv = f"{camera}:SPECT-COM"
|
||||
std_pv = f"{camera}:SPECT-RMS"
|
||||
skew_pv = f"{camera}:SPECT-SKEW"
|
||||
iqr_pv = f"{camera}:SPECT-IQR"
|
||||
res_pv = f"{camera}:SPECT-RES"
|
||||
|
||||
# Base PVs for update thread (order matters)
|
||||
base_pv_names = [
|
||||
e_int, center_pv, fwhm_pv, fit_rms_pv,
|
||||
fit_res_pv, com_pv, std_pv, skew_pv, iqr_pv, res_pv
|
||||
fit_res_pv, fit_spec_pv, com_pv, std_pv, skew_pv, iqr_pv, res_pv
|
||||
]
|
||||
|
||||
# Running-average configuration
|
||||
global_ravg_length = params.get('RAVG_length', global_ravg_length)
|
||||
# Build list of running-average PVs (exclude e_int, e_axis, processing_parameters)
|
||||
exclude = {
|
||||
e_int,
|
||||
e_axis,
|
||||
f"{camera}:processing_parameters"
|
||||
}
|
||||
exclude = {e_int, e_axis, f"{camera}:processing_parameters"}
|
||||
ravg_base = [pv for pv in base_pv_names if pv not in exclude]
|
||||
ravg_pv_names = [pv + '-RAVG' for pv in ravg_base]
|
||||
|
||||
# All PVs (original + running average)
|
||||
all_pv_names = base_pv_names + ravg_pv_names
|
||||
# N-shot average configuration
|
||||
global_avg_length = params.get('avg_nshots', global_avg_length)
|
||||
avg_buffer = deque(maxlen=global_avg_length)
|
||||
# Define PVs for N-shot average statistics and spectra
|
||||
avg_pv_names = [
|
||||
f"{camera}:AVG-FIT-COM",
|
||||
f"{camera}:AVG-FIT-FWHM",
|
||||
f"{camera}:AVG-FIT-RMS",
|
||||
f"{camera}:AVG-FIT-RES",
|
||||
f"{camera}:AVG-SPECT-COM",
|
||||
f"{camera}:AVG-SPECT-RMS",
|
||||
f"{camera}:AVG-SPECT-SKEW",
|
||||
f"{camera}:AVG-SPECT-IQR",
|
||||
f"{camera}:AVG-SPECT-RES",
|
||||
f"{camera}:AVG-SPECTRUM_Y",
|
||||
f"{camera}:AVG-FIT-SPECTRUM_Y"
|
||||
]
|
||||
|
||||
# All PVs (original + running average + N-shot average)
|
||||
all_pv_names = base_pv_names + ravg_pv_names + avg_pv_names
|
||||
|
||||
# Start background thread for PV updates
|
||||
thread = Thread(target=update_PVs, args=(buffer, *all_pv_names), daemon=True)
|
||||
@@ -109,15 +130,18 @@ def initialize(params):
|
||||
def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata=None, background=None):
|
||||
"""
|
||||
Main entrypoint: subtract background, crop ROI, smooth, fit Gaussian,
|
||||
compute metrics, queue PV updates (with running averages for skew and IQR).
|
||||
Returns a dict of processed PV values (original channels only).
|
||||
compute metrics, N-shot average, queue PV updates (with running averages).
|
||||
Returns a dict of processed PV values (original and average channels).
|
||||
"""
|
||||
global initialized, sent_pid, channel_pv_names, global_ravg_length, ravg_buffers
|
||||
global initialized, sent_pid, channel_pv_names
|
||||
global global_ravg_length, ravg_buffers, avg_buffer
|
||||
try:
|
||||
if not initialized:
|
||||
initialize(parameters)
|
||||
initialized = True
|
||||
|
||||
camera = parameters["camera_name"]
|
||||
|
||||
# Dynamic ROI and axis PV read
|
||||
ymin_pv, ymax_pv, axis_pv = create_thread_pvs(channel_pv_names)
|
||||
if ymin_pv and ymin_pv.connected:
|
||||
@@ -128,11 +152,7 @@ def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata
|
||||
if not (axis_pv and axis_pv.connected):
|
||||
_logger.warning("Energy axis not connected")
|
||||
return None
|
||||
axis = axis_pv.value
|
||||
if len(axis) < image.shape[1]:
|
||||
_logger.warning("Energy axis length %d < image width %d", len(axis), image.shape[1])
|
||||
return None
|
||||
axis = axis[:image.shape[1]]
|
||||
axis = axis_pv.value[:image.shape[1]]
|
||||
|
||||
# Preprocess image
|
||||
proc_img = image.astype(np.float32) - np.float32(parameters.get("pixel_bkg", 0))
|
||||
@@ -140,51 +160,36 @@ def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata
|
||||
|
||||
# Background image
|
||||
bg_img = parameters.pop('background_data', None)
|
||||
if not (isinstance(bg_img, np.ndarray) and bg_img.shape == proc_img.shape):
|
||||
bg_img = None
|
||||
else:
|
||||
bg_img = bg_img.astype(np.float32)
|
||||
bg_img = bg_img.astype(np.float32) if isinstance(bg_img, np.ndarray) and bg_img.shape == proc_img.shape else None
|
||||
|
||||
# Crop ROI
|
||||
ymin, ymax = int(global_roi[0]), int(global_roi[1])
|
||||
ymin, ymax = map(int, global_roi)
|
||||
if 0 <= ymin < ymax <= nrows:
|
||||
proc_img = proc_img[ymin:ymax, :]
|
||||
if bg_img is not None:
|
||||
bg_img = bg_img[ymin:ymax, :]
|
||||
|
||||
# Extract spectrum
|
||||
# Extract spectrum and fitted spectrum
|
||||
spectrum = get_spectrum(proc_img, bg_img) if bg_img is not None else np.sum(proc_img, axis=0)
|
||||
|
||||
# Smooth
|
||||
smoothed = scipy.signal.savgol_filter(spectrum, 51, 3)
|
||||
|
||||
# Noise check and fit Gaussian
|
||||
minimum, maximum = smoothed.min(), smoothed.max()
|
||||
amplitude = maximum - minimum
|
||||
skip = amplitude <= nrows * 1.5
|
||||
offset, amp_fit, center, sigma = functions.gauss_fit_psss(
|
||||
smoothed[::2], axis[::2], offset=minimum,
|
||||
amplitude=amplitude, skip=skip, maxfev=10
|
||||
smoothed[::2], axis[::2], offset=minimum, amplitude=amplitude, skip=skip, maxfev=10
|
||||
)
|
||||
# Reconstruct fitted curve
|
||||
fit_spectrum = offset + amp_fit * np.exp(-((axis - center)**2) / (2 * sigma**2))
|
||||
|
||||
# Compute normalized spectrum weights
|
||||
# Moments
|
||||
sm_norm = smoothed / np.sum(smoothed)
|
||||
|
||||
# Statistical moments
|
||||
spect_com = np.sum(axis * sm_norm)
|
||||
spect_std = np.sqrt(np.sum((axis - spect_com)**2 * sm_norm))
|
||||
spect_com = np.sum(axis * sm_norm)
|
||||
spect_std = np.sqrt(np.sum((axis - spect_com)**2 * sm_norm))
|
||||
spect_skew = np.sum((axis - spect_com)**3 * sm_norm) / (spect_std**3)
|
||||
cum = np.cumsum(sm_norm); e25 = np.interp(0.25, cum, axis); e75 = np.interp(0.75, cum, axis)
|
||||
spect_iqr = e75 - e25; spect_sum = np.sum(spectrum)
|
||||
|
||||
# Interquartile width (IQR)
|
||||
cum = np.cumsum(sm_norm)
|
||||
e25 = np.interp(0.25, cum, axis)
|
||||
e75 = np.interp(0.75, cum, axis)
|
||||
spect_iqr = e75 - e25
|
||||
|
||||
spect_sum = np.sum(spectrum)
|
||||
|
||||
camera = parameters["camera_name"]
|
||||
# Original result dict
|
||||
# Original result
|
||||
result = {
|
||||
parameters["e_int_name"]: spectrum,
|
||||
parameters["e_axis_name"]: axis,
|
||||
@@ -193,31 +198,59 @@ def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata
|
||||
f"{camera}:FIT-FWHM": np.float64(2.355 * sigma),
|
||||
f"{camera}:FIT-RMS": np.float64(sigma),
|
||||
f"{camera}:FIT-RES": np.float64(2.355 * sigma / center * 1000),
|
||||
f"{camera}:FIT-SPECTRUM_Y": fit_spectrum,
|
||||
f"{camera}:SPECT-COM": spect_com,
|
||||
f"{camera}:SPECT-RMS": spect_std,
|
||||
f"{camera}:SPECT-SKEW": spect_skew,
|
||||
f"{camera}:SPECT-IQR": spect_iqr,
|
||||
# Use IQR for relative spread instead of std
|
||||
f"{camera}:SPECT-RES": np.float64(spect_iqr / spect_com * 1000),
|
||||
f"{camera}:processing_parameters": json.dumps({"roi": global_roi})
|
||||
}
|
||||
|
||||
# Prepare full values for PV update (including running averages)
|
||||
exclude = {
|
||||
parameters["e_int_name"],
|
||||
parameters["e_axis_name"],
|
||||
f"{camera}:processing_parameters"
|
||||
}
|
||||
# Rolling averages
|
||||
exclude = {parameters["e_int_name"], parameters["e_axis_name"], f"{camera}:processing_parameters"}
|
||||
ravg_results = {}
|
||||
for base_pv in (pv for pv in base_pv_names if pv not in exclude):
|
||||
buf = ravg_buffers.setdefault(base_pv, deque(maxlen=global_ravg_length))
|
||||
buf.append(result.get(base_pv))
|
||||
buf.append(result[base_pv])
|
||||
ravg_results[f"{base_pv}-RAVG"] = np.mean(buf)
|
||||
|
||||
# Merge for PV write
|
||||
full_results = {**result, **ravg_results}
|
||||
# N-shot average and average fitted spectrum
|
||||
avg_buffer.append(spectrum)
|
||||
avg_spectrum = np.mean(np.stack(avg_buffer), axis=0)
|
||||
fit_avg = offset + amp_fit * np.exp(-((axis - center)**2) / (2 * sigma**2)) # using avg fit params below
|
||||
sm_avg = scipy.signal.savgol_filter(avg_spectrum, 51, 3)
|
||||
min_a, max_a = sm_avg.min(), sm_avg.max()
|
||||
amp_a = max_a - min_a; skip_a = amp_a <= nrows * 1.5
|
||||
offs_a, amp_fit_a, center_a, sigma_a = functions.gauss_fit_psss(
|
||||
sm_avg[::2], axis[::2], offset=min_a, amplitude=amp_a, skip=skip_a, maxfev=10
|
||||
)
|
||||
fit_avg_spectrum = offs_a + amp_fit_a * np.exp(-((axis - center_a)**2) / (2 * sigma_a**2))
|
||||
# Average moments
|
||||
sm_norm_a = sm_avg / np.sum(sm_avg)
|
||||
spect_com_a = np.sum(axis * sm_norm_a)
|
||||
spect_std_a = np.sqrt(np.sum((axis - spect_com_a)**2 * sm_norm_a))
|
||||
spect_skew_a = np.sum((axis - spect_com_a)**3 * sm_norm_a) / (spect_std_a**3)
|
||||
cum_a = np.cumsum(sm_norm_a); e25_a = np.interp(0.25, cum_a, axis); e75_a = np.interp(0.75, cum_a, axis)
|
||||
spect_iqr_a = e75_a - e25_a
|
||||
spect_res_a = spect_iqr_a / spect_com_a * 1000
|
||||
|
||||
# Queue PV update if new pulse
|
||||
avg_results = {
|
||||
f"{camera}:AVG-FIT-COM": np.float64(center_a),
|
||||
f"{camera}:AVG-FIT-FWHM": np.float64(2.355 * sigma_a),
|
||||
f"{camera}:AVG-FIT-RMS": np.float64(sigma_a),
|
||||
f"{camera}:AVG-FIT-RES": np.float64(2.355 * sigma_a / center_a * 1000),
|
||||
f"{camera}:AVG-SPECT-COM": spect_com_a,
|
||||
f"{camera}:AVG-SPECT-RMS": spect_std_a,
|
||||
f"{camera}:AVG-SPECT-SKEW": spect_skew_a,
|
||||
f"{camera}:AVG-SPECT-IQR": spect_iqr_a,
|
||||
f"{camera}:AVG-SPECT-RES": np.float64(spect_res_a),
|
||||
f"{camera}:AVG-SPECTRUM_Y": avg_spectrum,
|
||||
f"{camera}:AVG-FIT-SPECTRUM_Y": fit_avg_spectrum
|
||||
}
|
||||
|
||||
# Merge and queue
|
||||
full_results = {**result, **ravg_results, **avg_results}
|
||||
if epics_lock.acquire(False):
|
||||
try:
|
||||
if pulse_id > sent_pid:
|
||||
|
||||
@@ -0,0 +1,234 @@
|
||||
from logging import getLogger
|
||||
from cam_server.pipeline.data_processing import functions
|
||||
from cam_server.utils import create_thread_pvs, epics_lock
|
||||
from collections import deque
|
||||
import json
|
||||
import numpy as np
|
||||
import scipy.signal
|
||||
import numba
|
||||
import time
|
||||
import sys
|
||||
from threading import Thread
|
||||
|
||||
# Configure Numba to use multiple threads
|
||||
numba.set_num_threads(4)
|
||||
|
||||
_logger = getLogger(__name__)
|
||||
|
||||
# Shared state globals
|
||||
global_roi = [0, 0]
|
||||
initialized = False
|
||||
sent_pid = -1
|
||||
buffer = deque(maxlen=5)
|
||||
channel_pv_names = None
|
||||
base_pv_names = []
|
||||
all_pv_names = []
|
||||
global_ravg_length = 100
|
||||
ravg_buffers = {}
|
||||
|
||||
@numba.njit(parallel=False)
|
||||
def get_spectrum(image, background):
|
||||
"""Compute background-subtracted spectrum via row-wise summation."""
|
||||
y, x = image.shape
|
||||
profile = np.zeros(x, dtype=np.float64)
|
||||
for i in numba.prange(y):
|
||||
for j in range(x):
|
||||
profile[j] += image[i, j] - background[i, j]
|
||||
return profile
|
||||
|
||||
|
||||
def update_PVs(buffer, *pv_names):
|
||||
"""Continuously read from buffer and write to EPICS PVs."""
|
||||
pvs = create_thread_pvs(list(pv_names))
|
||||
while True:
|
||||
time.sleep(0.1)
|
||||
try:
|
||||
rec = buffer.popleft()
|
||||
except IndexError:
|
||||
continue
|
||||
try:
|
||||
for pv, val in zip(pvs, rec):
|
||||
if pv and pv.connected and (val is not None):
|
||||
pv.put(val)
|
||||
except Exception:
|
||||
_logger.exception("Error updating channels")
|
||||
|
||||
|
||||
def initialize(params):
|
||||
"""Initialize PV names, running-average settings, and launch update thread."""
|
||||
global channel_pv_names, base_pv_names, all_pv_names, global_ravg_length
|
||||
|
||||
camera = params["camera_name"]
|
||||
e_int = params["e_int_name"]
|
||||
e_axis = params["e_axis_name"]
|
||||
|
||||
# Fit/result PV names
|
||||
center_pv = f"{camera}:FIT-COM"
|
||||
fwhm_pv = f"{camera}:FIT-FWHM"
|
||||
fit_rms_pv = f"{camera}:FIT-RMS"
|
||||
fit_res_pv = f"{camera}:FIT-RES"
|
||||
|
||||
# ROI PVs for dynamic read
|
||||
ymin_pv = f"{camera}:SPC_ROI_YMIN"
|
||||
ymax_pv = f"{camera}:SPC_ROI_YMAX"
|
||||
axis_pv = e_axis
|
||||
channel_pv_names = [ymin_pv, ymax_pv, axis_pv]
|
||||
|
||||
# Spectrum statistical PV names
|
||||
com_pv = f"{camera}:SPECT-COM"
|
||||
std_pv = f"{camera}:SPECT-RMS"
|
||||
skew_pv = f"{camera}:SPECT-SKEW"
|
||||
iqr_pv = f"{camera}:SPECT-IQR"
|
||||
res_pv = f"{camera}:SPECT-RES" # will use IQR-based calc
|
||||
|
||||
# Base PVs for update thread (order matters)
|
||||
base_pv_names = [
|
||||
e_int, center_pv, fwhm_pv, fit_rms_pv,
|
||||
fit_res_pv, com_pv, std_pv, skew_pv, iqr_pv, res_pv
|
||||
]
|
||||
|
||||
# Running-average configuration
|
||||
global_ravg_length = params.get('RAVG_length', global_ravg_length)
|
||||
# Build list of running-average PVs (exclude e_int, e_axis, processing_parameters)
|
||||
exclude = {
|
||||
e_int,
|
||||
e_axis,
|
||||
f"{camera}:processing_parameters"
|
||||
}
|
||||
ravg_base = [pv for pv in base_pv_names if pv not in exclude]
|
||||
ravg_pv_names = [pv + '-RAVG' for pv in ravg_base]
|
||||
|
||||
# All PVs (original + running average)
|
||||
all_pv_names = base_pv_names + ravg_pv_names
|
||||
|
||||
# Start background thread for PV updates
|
||||
thread = Thread(target=update_PVs, args=(buffer, *all_pv_names), daemon=True)
|
||||
thread.start()
|
||||
|
||||
|
||||
def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata=None, background=None):
|
||||
"""
|
||||
Main entrypoint: subtract background, crop ROI, smooth, fit Gaussian,
|
||||
compute metrics, queue PV updates (with running averages for skew and IQR).
|
||||
Returns a dict of processed PV values (original channels only).
|
||||
"""
|
||||
global initialized, sent_pid, channel_pv_names, global_ravg_length, ravg_buffers
|
||||
try:
|
||||
if not initialized:
|
||||
initialize(parameters)
|
||||
initialized = True
|
||||
|
||||
# Dynamic ROI and axis PV read
|
||||
ymin_pv, ymax_pv, axis_pv = create_thread_pvs(channel_pv_names)
|
||||
if ymin_pv and ymin_pv.connected:
|
||||
global_roi[0] = ymin_pv.value
|
||||
if ymax_pv and ymax_pv.connected:
|
||||
global_roi[1] = ymax_pv.value
|
||||
|
||||
if not (axis_pv and axis_pv.connected):
|
||||
_logger.warning("Energy axis not connected")
|
||||
return None
|
||||
axis = axis_pv.value
|
||||
if len(axis) < image.shape[1]:
|
||||
_logger.warning("Energy axis length %d < image width %d", len(axis), image.shape[1])
|
||||
return None
|
||||
axis = axis[:image.shape[1]]
|
||||
|
||||
# Preprocess image
|
||||
proc_img = image.astype(np.float32) - np.float32(parameters.get("pixel_bkg", 0))
|
||||
nrows, _ = proc_img.shape
|
||||
|
||||
# Background image
|
||||
bg_img = parameters.pop('background_data', None)
|
||||
if not (isinstance(bg_img, np.ndarray) and bg_img.shape == proc_img.shape):
|
||||
bg_img = None
|
||||
else:
|
||||
bg_img = bg_img.astype(np.float32)
|
||||
|
||||
# Crop ROI
|
||||
ymin, ymax = int(global_roi[0]), int(global_roi[1])
|
||||
if 0 <= ymin < ymax <= nrows:
|
||||
proc_img = proc_img[ymin:ymax, :]
|
||||
if bg_img is not None:
|
||||
bg_img = bg_img[ymin:ymax, :]
|
||||
|
||||
# Extract spectrum
|
||||
spectrum = get_spectrum(proc_img, bg_img) if bg_img is not None else np.sum(proc_img, axis=0)
|
||||
|
||||
# Smooth
|
||||
smoothed = scipy.signal.savgol_filter(spectrum, 51, 3)
|
||||
|
||||
# Noise check and fit Gaussian
|
||||
minimum, maximum = smoothed.min(), smoothed.max()
|
||||
amplitude = maximum - minimum
|
||||
skip = amplitude <= nrows * 1.5
|
||||
offset, amp_fit, center, sigma = functions.gauss_fit_psss(
|
||||
smoothed[::2], axis[::2], offset=minimum,
|
||||
amplitude=amplitude, skip=skip, maxfev=10
|
||||
)
|
||||
|
||||
# Compute normalized spectrum weights
|
||||
sm_norm = smoothed / np.sum(smoothed)
|
||||
|
||||
# Statistical moments
|
||||
spect_com = np.sum(axis * sm_norm)
|
||||
spect_std = np.sqrt(np.sum((axis - spect_com)**2 * sm_norm))
|
||||
spect_skew = np.sum((axis - spect_com)**3 * sm_norm) / (spect_std**3)
|
||||
|
||||
# Interquartile width (IQR)
|
||||
cum = np.cumsum(sm_norm)
|
||||
e25 = np.interp(0.25, cum, axis)
|
||||
e75 = np.interp(0.75, cum, axis)
|
||||
spect_iqr = e75 - e25
|
||||
|
||||
spect_sum = np.sum(spectrum)
|
||||
|
||||
camera = parameters["camera_name"]
|
||||
# Original result dict
|
||||
result = {
|
||||
parameters["e_int_name"]: spectrum,
|
||||
parameters["e_axis_name"]: axis,
|
||||
f"{camera}:SPECTRUM_Y_SUM": spect_sum,
|
||||
f"{camera}:FIT-COM": np.float64(center),
|
||||
f"{camera}:FIT-FWHM": np.float64(2.355 * sigma),
|
||||
f"{camera}:FIT-RMS": np.float64(sigma),
|
||||
f"{camera}:FIT-RES": np.float64(2.355 * sigma / center * 1000),
|
||||
f"{camera}:SPECT-COM": spect_com,
|
||||
f"{camera}:SPECT-RMS": spect_std,
|
||||
f"{camera}:SPECT-SKEW": spect_skew,
|
||||
f"{camera}:SPECT-IQR": spect_iqr,
|
||||
# Use IQR for relative spread instead of std
|
||||
f"{camera}:SPECT-RES": np.float64(spect_iqr / spect_com * 1000),
|
||||
f"{camera}:processing_parameters": json.dumps({"roi": global_roi})
|
||||
}
|
||||
|
||||
# Prepare full values for PV update (including running averages)
|
||||
exclude = {
|
||||
parameters["e_int_name"],
|
||||
parameters["e_axis_name"],
|
||||
f"{camera}:processing_parameters"
|
||||
}
|
||||
ravg_results = {}
|
||||
for base_pv in (pv for pv in base_pv_names if pv not in exclude):
|
||||
buf = ravg_buffers.setdefault(base_pv, deque(maxlen=global_ravg_length))
|
||||
buf.append(result.get(base_pv))
|
||||
ravg_results[f"{base_pv}-RAVG"] = np.mean(buf)
|
||||
|
||||
# Merge for PV write
|
||||
full_results = {**result, **ravg_results}
|
||||
|
||||
# Queue PV update if new pulse
|
||||
if epics_lock.acquire(False):
|
||||
try:
|
||||
if pulse_id > sent_pid:
|
||||
sent_pid = pulse_id
|
||||
entry = tuple(full_results.get(pv) for pv in all_pv_names)
|
||||
buffer.append(entry)
|
||||
finally:
|
||||
epics_lock.release()
|
||||
|
||||
return full_results
|
||||
|
||||
except Exception as ex:
|
||||
_logger.warning("Processing error: %s", ex)
|
||||
return {}
|
||||
@@ -0,0 +1,283 @@
|
||||
from logging import getLogger
|
||||
from cam_server.pipeline.data_processing import functions
|
||||
from cam_server.utils import epics_lock, create_thread_pvs
|
||||
from collections import deque
|
||||
import threading
|
||||
import json
|
||||
import numpy as np
|
||||
import scipy.signal
|
||||
import numba
|
||||
import time
|
||||
from threading import Thread
|
||||
|
||||
# Configure Numba to use multiple threads and parallelize
|
||||
numba.set_num_threads(4)
|
||||
|
||||
_logger = getLogger(__name__)
|
||||
|
||||
# Shared state globals
|
||||
global_roi = [0, 0]
|
||||
last_roi = None
|
||||
params_json = None
|
||||
initialized = False
|
||||
sent_pid = -1
|
||||
buffer = deque(maxlen=5)
|
||||
channel_pv_names = []
|
||||
base_pv_names = []
|
||||
all_pv_names = []
|
||||
global_ravg_length = 100
|
||||
# Running-average state
|
||||
global ravg_buffers, ravg_sum_map, ravg_lock
|
||||
ravg_buffers = {}
|
||||
ravg_sum_map = {}
|
||||
ravg_lock = threading.Lock()
|
||||
# N-shot average state
|
||||
global global_avg_n, avg_buffer, avg_sum
|
||||
global_avg_n = 100
|
||||
avg_buffer = deque()
|
||||
avg_sum = None
|
||||
# Cached PV handles
|
||||
pv_handles = None
|
||||
|
||||
# Update thread function
|
||||
def update_PVs():
|
||||
"""Continuously pop and write to EPICS PVs using cached handles."""
|
||||
global pv_handles, buffer
|
||||
while True:
|
||||
try:
|
||||
entry = buffer.popleft()
|
||||
except IndexError:
|
||||
time.sleep(0.001)
|
||||
continue
|
||||
for pv, val in zip(pv_handles, entry):
|
||||
if pv and pv.connected and val is not None:
|
||||
pv.put(val)
|
||||
|
||||
|
||||
def initialize(params):
|
||||
"""Initialize PV names, caches, and start update thread."""
|
||||
global channel_pv_names, base_pv_names, all_pv_names
|
||||
global global_ravg_length, global_avg_n, avg_buffer, avg_sum, pv_handles
|
||||
|
||||
camera = params["camera_name"]
|
||||
e_int = params["e_int_name"]
|
||||
e_axis = params["e_axis_name"]
|
||||
|
||||
# Fit/result PV names
|
||||
center_pv = f"{camera}:FIT-COM"
|
||||
fwhm_pv = f"{camera}:FIT-FWHM"
|
||||
fit_rms_pv = f"{camera}:FIT-RMS"
|
||||
fit_res_pv = f"{camera}:FIT-RES"
|
||||
fit_spec_pv = f"{camera}:FIT-SPECTRUM_Y"
|
||||
|
||||
# ROI PVs
|
||||
ymin_pv = f"{camera}:SPC_ROI_YMIN"
|
||||
ymax_pv = f"{camera}:SPC_ROI_YMAX"
|
||||
channel_pv_names = [ymin_pv, ymax_pv, e_axis]
|
||||
|
||||
# Stats PVs
|
||||
com_pv = f"{camera}:SPECT-COM"
|
||||
std_pv = f"{camera}:SPECT-RMS"
|
||||
skew_pv = f"{camera}:SPECT-SKEW"
|
||||
iqr_pv = f"{camera}:SPECT-IQR"
|
||||
res_pv = f"{camera}:SPECT-RES"
|
||||
|
||||
# Base PV list
|
||||
base_pv_names = [
|
||||
e_int, center_pv, fwhm_pv, fit_rms_pv,
|
||||
fit_res_pv, fit_spec_pv, com_pv, std_pv, skew_pv, iqr_pv, res_pv
|
||||
]
|
||||
|
||||
# Running-average PV names
|
||||
global_ravg_length = params.get('RAVG_length', global_ravg_length)
|
||||
exclude = {e_int, e_axis, f"{camera}:processing_parameters"}
|
||||
ravg_base = [pv for pv in base_pv_names if pv not in exclude]
|
||||
ravg_pv_names = [pv + '-RAVG' for pv in ravg_base]
|
||||
|
||||
# N-shot average settings
|
||||
global_avg_n = params.get('avg_nshots', global_avg_n)
|
||||
avg_buffer = deque(maxlen=global_avg_n)
|
||||
avg_sum = None
|
||||
avg_pv_names = [
|
||||
f"{camera}:AVG-FIT-COM",
|
||||
f"{camera}:AVG-FIT-FWHM",
|
||||
f"{camera}:AVG-FIT-RMS",
|
||||
f"{camera}:AVG-FIT-RES",
|
||||
f"{camera}:AVG-SPECT-COM",
|
||||
f"{camera}:AVG-SPECT-RMS",
|
||||
f"{camera}:AVG-SPECT-SKEW",
|
||||
f"{camera}:AVG-SPECT-IQR",
|
||||
f"{camera}:AVG-SPECT-RES",
|
||||
f"{camera}:AVG-SPECTRUM_Y",
|
||||
f"{camera}:AVG-FIT-SPECTRUM_Y"
|
||||
]
|
||||
|
||||
# All PVs
|
||||
all_pv_names = base_pv_names + ravg_pv_names + avg_pv_names + [f"{camera}:processing_parameters"]
|
||||
|
||||
# Cache PV handles
|
||||
pv_handles = create_thread_pvs(all_pv_names)
|
||||
|
||||
# Start update thread
|
||||
thread = Thread(target=update_PVs, daemon=True)
|
||||
thread.start()
|
||||
|
||||
|
||||
def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata=None, background=None):
|
||||
"""
|
||||
Fast processing: vectorized spectrum, incremental averages, cached PVs.
|
||||
"""
|
||||
global initialized, sent_pid, buffer
|
||||
global global_roi, last_roi, params_json
|
||||
global ravg_buffers, ravg_sum_map, ravg_lock
|
||||
global avg_buffer, avg_sum, global_avg_n, all_pv_names
|
||||
|
||||
if not initialized:
|
||||
initialize(parameters)
|
||||
initialized = True
|
||||
|
||||
camera = parameters["camera_name"]
|
||||
|
||||
# Read ROI once and detect change
|
||||
ymin, ymax = global_roi
|
||||
new_ymin = parameters.get('roi_ymin', ymin)
|
||||
new_ymax = parameters.get('roi_ymax', ymax)
|
||||
if (new_ymin, new_ymax) != (ymin, ymax):
|
||||
global_roi[:] = [new_ymin, new_ymax]
|
||||
roi_changed = True
|
||||
else:
|
||||
roi_changed = False
|
||||
|
||||
# Preprocess image
|
||||
img = image.astype(np.float32)
|
||||
pixel_bkg = parameters.get("pixel_bkg", 0)
|
||||
proc_img = img - pixel_bkg
|
||||
nrows, _ = proc_img.shape
|
||||
|
||||
# Background subtraction
|
||||
bg = parameters.pop('background_data', None)
|
||||
if isinstance(bg, np.ndarray) and bg.shape == proc_img.shape:
|
||||
proc_img -= bg.astype(np.float32)
|
||||
|
||||
# Crop ROI
|
||||
ymin_i, ymax_i = int(global_roi[0]), int(global_roi[1])
|
||||
if 0 <= ymin_i < ymax_i <= nrows:
|
||||
proc_img = proc_img[ymin_i:ymax_i, :]
|
||||
|
||||
# Spectrum via vector sum
|
||||
spectrum = np.sum(proc_img, axis=0)
|
||||
# Smooth
|
||||
smoothed = scipy.signal.savgol_filter(spectrum, 51, 3)
|
||||
# Fit Gaussian
|
||||
minimum, maximum = smoothed.min(), smoothed.max()
|
||||
amplitude = maximum - minimum
|
||||
skip = amplitude <= nrows * 1.5
|
||||
offset, amp_fit, center, sigma = functions.gauss_fit_psss(
|
||||
smoothed[::2], x_axis[:len(smoothed)][::2], offset=minimum,
|
||||
amplitude=amplitude, skip=skip, maxfev=10
|
||||
)
|
||||
fit_spectrum = offset + amp_fit * np.exp(-((x_axis[:len(smoothed)] - center)**2) / (2 * sigma**2))
|
||||
|
||||
# Compute stats
|
||||
sm_norm = smoothed / np.sum(smoothed)
|
||||
spect_com = np.dot(x_axis[:len(sm_norm)], sm_norm)
|
||||
spect_std = np.sqrt(np.dot((x_axis[:len(sm_norm)] - spect_com)**2, sm_norm))
|
||||
spect_skew = np.dot((x_axis[:len(sm_norm)] - spect_com)**3, sm_norm) / (spect_std**3)
|
||||
cum = np.cumsum(sm_norm)
|
||||
e25 = np.interp(0.25, cum, x_axis[:len(cum)])
|
||||
e75 = np.interp(0.75, cum, x_axis[:len(cum)])
|
||||
spect_iqr = e75 - e25
|
||||
spect_sum = spectrum.sum()
|
||||
|
||||
# Original result
|
||||
result = {
|
||||
parameters["e_int_name"]: spectrum,
|
||||
parameters["e_axis_name"]: x_axis[:len(spectrum)],
|
||||
f"{camera}:SPECTRUM_Y_SUM": spect_sum,
|
||||
f"{camera}:FIT-COM": np.float64(center),
|
||||
f"{camera}:FIT-FWHM": np.float64(2.355 * sigma),
|
||||
f"{camera}:FIT-RMS": np.float64(sigma),
|
||||
f"{camera}:FIT-RES": np.float64(2.355 * sigma / center * 1000),
|
||||
f"{camera}:FIT-SPECTRUM_Y": fit_spectrum,
|
||||
f"{camera}:SPECT-COM": spect_com,
|
||||
f"{camera}:SPECT-RMS": spect_std,
|
||||
f"{camera}:SPECT-SKEW": spect_skew,
|
||||
f"{camera}:SPECT-IQR": spect_iqr,
|
||||
f"{camera}:SPECT-RES": np.float64(spect_iqr / spect_com * 1000)
|
||||
}
|
||||
# Update JSON PV only on ROI change
|
||||
if roi_changed:
|
||||
params_json = json.dumps({"roi": global_roi})
|
||||
result[f"{camera}:processing_parameters"] = params_json
|
||||
|
||||
# Running averages (thread-safe, incremental sum)
|
||||
ravg_results = {}
|
||||
exclude = {parameters["e_int_name"], parameters["e_axis_name"], f"{camera}:processing_parameters"}
|
||||
with ravg_lock:
|
||||
for pv in base_pv_names:
|
||||
if pv not in exclude:
|
||||
buf = ravg_buffers.setdefault(pv, deque(maxlen=global_ravg_length))
|
||||
sum_val = ravg_sum_map.get(pv, 0.0)
|
||||
if len(buf) == buf.maxlen:
|
||||
old = buf.popleft()
|
||||
sum_val -= old
|
||||
buf.append(result[pv])
|
||||
sum_val += result[pv]
|
||||
ravg_sum_map[pv] = sum_val
|
||||
ravg_results[pv + "-RAVG"] = sum_val / len(buf)
|
||||
|
||||
# N-shot average (incremental, no full-stack)
|
||||
global avg_sum
|
||||
if avg_sum is None:
|
||||
avg_sum = np.zeros_like(spectrum)
|
||||
if len(avg_buffer) == global_avg_n:
|
||||
old = avg_buffer.popleft()
|
||||
avg_sum -= old
|
||||
avg_buffer.append(spectrum)
|
||||
avg_sum += spectrum
|
||||
avg_spectrum = avg_sum / len(avg_buffer)
|
||||
# Fit and stats on avg_spectrum
|
||||
sm_avg = scipy.signal.savgol_filter(avg_spectrum, 51, 3)
|
||||
min_a, max_a = sm_avg.min(), sm_avg.max()
|
||||
amp_a = max_a - min_a
|
||||
skip_a = amp_a <= nrows * 1.5
|
||||
offs_a, amp_fit_a, center_a, sigma_a = functions.gauss_fit_psss(
|
||||
sm_avg[::2], x_axis[:len(sm_avg)][::2], offset=min_a,
|
||||
amplitude=amp_a, skip=skip_a, maxfev=10
|
||||
)
|
||||
fit_avg_spectrum = offs_a + amp_fit_a * np.exp(-((x_axis[:len(sm_avg)] - center_a)**2) / (2 * sigma_a**2))
|
||||
sm_norm_a = sm_avg / np.sum(sm_avg)
|
||||
spect_com_a = np.dot(x_axis[:len(sm_norm_a)], sm_norm_a)
|
||||
spect_std_a = np.sqrt(np.dot((x_axis[:len(sm_norm_a)] - spect_com_a)**2, sm_norm_a))
|
||||
spect_skew_a = np.dot((x_axis[:len(sm_norm_a)] - spect_com_a)**3, sm_norm_a) / (spect_std_a**3)
|
||||
cum_a = np.cumsum(sm_norm_a)
|
||||
e25_a = np.interp(0.25, cum_a, x_axis[:len(cum_a)])
|
||||
e75_a = np.interp(0.75, cum_a, x_axis[:len(cum_a)])
|
||||
spect_iqr_a = e75_a - e25_a
|
||||
spect_res_a = spect_iqr_a / spect_com_a * 1000
|
||||
|
||||
avg_results = {
|
||||
f"{camera}:AVG-FIT-COM": np.float64(center_a),
|
||||
f"{camera}:AVG-FIT-FWHM": np.float64(2.355 * sigma_a),
|
||||
f"{camera}:AVG-FIT-RMS": np.float64(sigma_a),
|
||||
f"{camera}:AVG-FIT-RES": np.float64(2.355 * sigma_a / center_a * 1000),
|
||||
f"{camera}:AVG-SPECT-COM": spect_com_a,
|
||||
f"{camera}:AVG-SPECT-RMS": spect_std_a,
|
||||
f"{camera}:AVG-SPECT-SKEW": spect_skew_a,
|
||||
f"{camera}:AVG-SPECT-IQR": spect_iqr_a,
|
||||
f"{camera}:AVG-SPECT-RES": np.float64(spect_res_a),
|
||||
f"{camera}:AVG-SPECTRUM_Y": avg_spectrum,
|
||||
f"{camera}:AVG-FIT-SPECTRUM_Y": fit_avg_spectrum
|
||||
}
|
||||
|
||||
# Merge and queue for PV update
|
||||
full = {**result, **ravg_results, **avg_results}
|
||||
if epics_lock.acquire(False):
|
||||
try:
|
||||
if pulse_id > sent_pid:
|
||||
sent_pid = pulse_id
|
||||
buffer.append(tuple(full[pv] for pv in all_pv_names))
|
||||
finally:
|
||||
epics_lock.release()
|
||||
|
||||
return full
|
||||
Reference in New Issue
Block a user