From c948942a12de7134b7170b1f6943c94f1e472d7e Mon Sep 17 00:00:00 2001 From: gac-alvra Date: Thu, 24 Mar 2022 13:33:11 +0100 Subject: [PATCH] renamed tt.py -> SARES11-SPEC125-M1-tt.py; updates --- procprof/{tt.py => SARES11-SPEC125-M1-tt.py} | 179 +++++++++---------- 1 file changed, 89 insertions(+), 90 deletions(-) rename procprof/{tt.py => SARES11-SPEC125-M1-tt.py} (67%) diff --git a/procprof/tt.py b/procprof/SARES11-SPEC125-M1-tt.py similarity index 67% rename from procprof/tt.py rename to procprof/SARES11-SPEC125-M1-tt.py index 08e588a..77b25d8 100644 --- a/procprof/tt.py +++ b/procprof/SARES11-SPEC125-M1-tt.py @@ -7,9 +7,82 @@ from scipy.ndimage import gaussian_filter1d from scipy.ndimage import convolve1d import numpy as np +# Alvra spectral encoder constants/waveforms +px2fs = 1.91 # calibration from ... +lambdas = 467.55 + 0.07219*np.arange(0,2048) # calibration from 23-9-2020 +nus = 299792458 / (lambdas * 10**-9) # frequency space, uneven +nus_new = np.linspace(nus[0], nus[-1], num=2048, endpoint=True) # frequency space, even +filters = { + 'YAG': np.concatenate((np.ones(50),signal.tukey(40)[20:40], np.zeros(1978), np.zeros(2048))), # fourier filter for YAGS + 'SiN': np.concatenate((signal.tukey(40)[20:40], np.zeros(2028), np.zeros(2048))), # fourier filter for 5um SiN + 'SiN2': np.concatenate((signal.tukey(32)[16:32], np.zeros(2032), np.zeros(2048))), # fourier filter for 2um SiN + 'babyYAG': np.concatenate((signal.tukey(40)[20:40], np.zeros(2028), np.zeros(2048))), # baby timetool YAG filter + 'babyYAG2': np.concatenate((np.ones(50),signal.tukey(40)[20:40], np.zeros(1978), np.zeros(2048))) # baby timetool YAG +} + +# Functions for image analysis and spectral encoding processing +def get_roi_projection(image, roi, axis): + x_start, x_stop, y_start, y_stop = roi + cropped = image[x_start:x_stop, y_start:y_stop] + project = cropped.mean(axis=axis) + return project + + +def interpolate(xold, xnew, vals): + ''' + Interpolate vals from xold to xnew + ''' + interp = interp1d(xold, vals, kind='cubic') + return interp(xnew) + + + +def fourier_filter(vals, filt): + ''' + Fourier transform, filter, inverse fourier transform, take the real part + ''' + vals = np.hstack((vals, np.zeros_like(vals))) # pad + transformed = np.fft.fft(vals) + filtered = transformed * filt + inverse = np.fft.ifft(filtered) + invreal = 2 * np.real(inverse) + return invreal + + + +def edge(filter_name, backgrounds, signals, background_from_fit, peakback): + ''' + returns: + edge positions determined from argmax of peak traces, converted to fs + edge amplitudes determined from the amax of the peak traces + the actual peak traces, which are the derivative of signal traces (below) + the raw signal traces, without any processing or smoothing except for the Fourier filter applied to remove the etalon + ''' + ffilter = filters[filter_name] + + # background normalization + sig_norm = np.nan_to_num(signals / backgrounds) / background_from_fit + # interpolate to get evenly sampled in frequency space + sig_interp = interpolate(nus, nus_new, sig_norm) + # Fourier filter + sig_filtered = fourier_filter(sig_interp, ffilter) + # interpolate to get unevenly sampled in frequency space (back to original wavelength space) + sig_uninterp = interpolate(nus_new, nus, sig_filtered[..., 0:2048]) + # peak via the derivative + sig_deriv = gaussian_filter1d(sig_uninterp, 50, order=1) + sig_deriv -= peakback + peak_pos = 1024 - np.argmax(sig_deriv, axis=-1) + peak_pos *= px2fs + peak_amp = np.amax(sig_deriv, axis=-1) + + return peak_pos, peak_amp, sig_deriv, sig_uninterp + def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata=None): + ''' + takes an image, processes the signal, sends out the raw data and processed results + ''' image = image.astype(int) camera_name = parameters["camera_name"] @@ -23,32 +96,30 @@ def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata project_axis = parameters.get("project_axis", 0) threshold = parameters.get("threshold") - # maintain the structure of processing_parameters background_shape = None # maintain the structure of res projection_signal = projection_background = None - try: if background is not None: - background_shape = background.shape - image -= background - - if threshold is not None: - image -= threshold - image[image < 0] = 0 +# background_shape = background.shape +# image -= background + projection_background = get_roi_projection(background, roi_signal, project_axis) +# +# if threshold is not None: +# image -= threshold +# image[image < 0] = 0 if roi_signal is not None: projection_signal = get_roi_projection(image, roi_signal, project_axis) +# +# if roi_background is not None: +# projection_background = get_roi_projection(image, roi_background, project_axis) - if roi_background is not None: - projection_background = get_roi_projection(image, roi_background, project_axis) - - peak2, sig6, sig5gaussO1 = edge("YAG", projection_background, projection_signal, 1, 0) -# print("peak2", peak2) + peak_pos, peak_amp, sig_deriv, sig_uninterp = edge("YAG", projection_background, projection_signal, 1, 0) except Exception as e: lineno = sys.exc_info()[2].tb_lineno @@ -58,7 +129,6 @@ def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata else: status = "OK" - processing_parameters = { "image_shape": image.shape, "background_shape": background_shape, @@ -80,82 +150,11 @@ def process_image(image, pulse_id, timestamp, x_axis, y_axis, parameters, bsdata camera_name + ".processing_parameters": processing_parameters, camera_name + ".projection_signal": projection_signal, - camera_name + ".projection_background": projection_background + camera_name + ".projection_background": projection_background, + camera_name + ".edge_position": peak_pos, + camera_name + ".edge_amplitude": peak_amp, + camera_name + ".edge_derivative": sig_deriv, + camera_name + ".edge_raw": sig_uninterp } return res - - -def get_roi_projection(image, roi, axis): - x_start, x_stop, y_start, y_stop = roi - cropped = image[x_start:x_stop, y_start:y_stop] - project = cropped.mean(axis=axis) - return project - - -lambdas = 467.55 + 0.07219*np.arange(0,2048) # calibration from 23-9-2020 -nus = 299792458 / (lambdas * 10**-9) # frequency space, uneven -nus_new = np.linspace(nus[0], nus[-1], num=2048, endpoint=True) # frequency space, even -derivFilter = np.concatenate((signal.tukey(2000)[0:500], np.ones(2048-1000), signal.tukey(2000)[1500:2000])) -Heaviside = np.concatenate((np.zeros(100), np.ones(100))) - - -filters = { - "YAG": np.concatenate((np.ones(50),signal.tukey(40)[20:40], np.zeros(1978), np.zeros(2048))), # fourier filter for YAGS - "SiN": np.concatenate((signal.tukey(40)[20:40], np.zeros(2028), np.zeros(2048))), # fourier filter for 5um SiN - "SiN2": np.concatenate((signal.tukey(32)[16:32], np.zeros(2032), np.zeros(2048))), # fourier filter for 2um SiN - "babyYAG": np.concatenate((signal.tukey(40)[20:40], np.zeros(2028), np.zeros(2048))), # baby timetool YAG filter - "babyYAG2": np.concatenate((np.ones(50),signal.tukey(40)[20:40], np.zeros(1978), np.zeros(2048))) # baby timetool YAG -} - -def edge(filter_name, backgrounds, signals, background_from_fit, peakback): - """ - returns: - edge positions determined from argmax of peak traces - signal traces, should show a change in transmission near px 1024 if set up correctly - peak traces, which are the derivative of signal traces - """ - ffilter = filters[filter_name] - - # background normalization - sig_norm = np.nan_to_num(signals / backgrounds) / background_from_fit - - # interpolate to get evenly sampled in frequency space - sig_interp = interpolate(nus, nus_new, sig_norm) - - # Fourier filter - sig_filtered = fourier_filter(sig_interp, ffilter) - - # interpolate to get unevenly sampled in frequency space (back to original wavelength space) - sig_uninterp = interpolate(nus_new, nus, sig_filtered[..., 0:2048]) - - # peak via the derivative - sig_deriv = gaussian_filter1d(sig_uninterp, 50, order=1) - sig_deriv -= peakback - peak_pos = np.argmax(sig_deriv, axis=-1) - - return peak_pos, sig_deriv, sig_uninterp - - - -def fourier_filter(vals, filt): - """ - Fourier transform, filter, inverse fourier transform, take the real part - """ - vals = np.hstack((vals, np.zeros_like(vals))) # pad - transformed = np.fft.fft(vals) - filtered = transformed * filt - inverse = np.fft.ifft(filtered) - invreal = 2 * np.real(inverse) - return invreal - - -def interpolate(xold, xnew, vals): - """ - Interpolate vals from xold to xnew - """ - interp = interp1d(xold, vals, kind='cubic') - return interp(xnew) - - -