renamed tt.py -> SARES11-SPEC125-M1-tt.py; updates
This commit is contained in:
gac-alvra
@ -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)
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user