from coseda.logging_utils import log_print
import h5py
import matplotlib.pyplot as plt
import h5py
import numpy as np
from scipy.signal import find_peaks
[docs]
def find_strokes_by_xpos(input_path, comparison_window = 10, ignore_window = 10):
with h5py.File(input_path, 'r') as file:
stagepos_x = file['entry/data/stagepos_x_refined'][:]
rightstrokes = []
leftstrokes = []
frame_numbers = np.arange(len(stagepos_x))
# Determine if the first significant change is a rightstroke or a leftstroke
initial_frames = stagepos_x[:comparison_window + 1]
looking_for_rightstroke = initial_frames[0] < max(initial_frames)
i = 0
while i < len(stagepos_x) - comparison_window:
current_frame = stagepos_x[i]
next_frames = stagepos_x[i+1:i+1+comparison_window]
if looking_for_rightstroke:
if all(current_frame < frame for frame in next_frames):
if not rightstrokes or (leftstrokes and i - leftstrokes[-1] > ignore_window):
rightstrokes.append(i)
looking_for_rightstroke = False
# Move the index forward by comparison window
i += comparison_window
else:
i += 1
else:
if all(current_frame > frame for frame in next_frames):
if not leftstrokes or (rightstrokes and i - rightstrokes[-1] > ignore_window):
leftstrokes.append(i)
looking_for_rightstroke = True
# Move the index forward by comparison window
i += comparison_window
else:
i += 1
"""
# Plotting
plt.figure(figsize=(14, 7))
plt.plot(frame_numbers, stagepos_x, label='Stage Position X')
# Plotting the beginning of rightstrokes
plt.scatter(rightstrokes, stagepos_x[rightstrokes], color='red', label='Beginning of right stroke')
# Plotting the beginning of leftstrokes
plt.scatter(leftstrokes, stagepos_x[leftstrokes], color='green', label='Beginning of left stroke')
plt.title('Stage Position X vs Frame Number with Rightstrokes and Leftstrokes')
plt.xlabel('Frame Number')
plt.ylabel('Stage Position X')
plt.legend()
plt.show()
log_print(f'X-method: {len(rightstrokes) + len(leftstrokes)}')
"""
return rightstrokes, leftstrokes
[docs]
def find_strokes_by_ypos(input_path, threshold_multiplier = 10):
with h5py.File(input_path, 'r') as file:
stagepos_y = file['entry/data/stagepos_y'][:]
line_break_indices = []
line_break_steps = [] # To store the step sizes at line breaks
refined_line_break_indices = [0]
refined_line_break_steps = [] # To store the step sizes at line breaks
window_size = 10 # Number of frames to consider for calculating standard deviation
ignore_window = 10 # Number of frames to ignore subsequent line breaks after a detected line break
# Calculate step sizes (difference between consecutive Y coordinates)
steps = np.abs(np.diff(stagepos_y))
for i in range(window_size, len(steps)):
# Calculate standard deviation of the last 'window_size' steps
window_std = np.std(steps[i-window_size:i])
# Check if the current step is significantly larger than recent variability
if steps[i] > threshold_multiplier * window_std:
# Check if the last line break was more than ignore_window frames ago
if not line_break_indices or i - line_break_indices[-1] > ignore_window:
line_break_indices.append(i + 1) # +1 as we are working with steps (diff)
line_break_steps.append(steps[i])
# Calculate and print the median step size of line breaks
median_step_size = np.median(line_break_steps)
#print("Median Step Size of Line Breaks:", median_step_size)
# Calculate median number of frames between each line break
frames_between_breaks = np.diff(line_break_indices)
median_frames_between_breaks = np.median(frames_between_breaks)
#print("Median Number of Frames Between Line Breaks:", median_frames_between_breaks)
for i in range(len(steps)):
if steps[i] > 0.5 * median_step_size:
# Using half of the median number of frames between line breaks from the previous run as new ignore window
if not refined_line_break_indices or i - refined_line_break_indices[-1] > 0.1 * median_frames_between_breaks:
refined_line_break_indices.append(i + 1) # +1 as we are working with steps (diff)
refined_line_break_steps.append(steps[i])
return line_break_indices, refined_line_break_indices, median_step_size
[docs]
def compare_strokes(leftstrokes, rightstrokes, refined_line_breaks):
unmatched_leftstrokes = []
unmatched_rightstrokes = []
# Create a combined list of all the strokes found by X method
all_strokes_x = sorted(leftstrokes + rightstrokes)
# For each stroke found by X method, check if there is a match in the Y method results
for stroke_x in all_strokes_x:
# Check if there's a matching frame in the Y method within a ±15 frame range
if not any(abs(stroke_x - stroke_y) <= 15 for stroke_y in refined_line_breaks):
# If there's no match, add to the corresponding unmatched list
if stroke_x in leftstrokes:
unmatched_leftstrokes.append(stroke_x)
else:
unmatched_rightstrokes.append(stroke_x)
return unmatched_leftstrokes, unmatched_rightstrokes
[docs]
def plot_stagepos_y_with_breaks(input_path):
# Plotting line breaks based on y position for debugging
rightstrokes, leftstrokes = find_strokes_by_xpos(input_path)
xmethod = len(rightstrokes) + len(leftstrokes)
log_print(f'x-method: {len(xmethod)}')
line_breaks, refined_line_breaks, median_step_size = find_strokes_by_ypos(input_path)
with h5py.File(input_path, 'r') as file:
stagepos_y = file['entry/data/stagepos_y'][:]
log_print(f'y-method: {len(refined_line_breaks)}')
plt.figure(figsize=(14, 7))
plt.plot(stagepos_y, label='Stage Position Y')
# Plot the initial line breaks based on Y coordinate analysis
#for lb in line_breaks:
# plt.axvline(x=lb, color='orange', linestyle=':', label='Initial Y Breaks' if lb == line_breaks[0] else "")
unmatched_leftstrokes, unmatched_rightstrokes = compare_strokes(leftstrokes, rightstrokes, refined_line_breaks)
log_print("Unmatched leftstrokes (X method):", unmatched_leftstrokes)
log_print("Unmatched rightstrokes (X method):", unmatched_rightstrokes)
# Plot the refined line breaks based on Y coordinate analysis
for lb in refined_line_breaks:
plt.axvline(x=lb, color='green', linestyle='--', label='Refined Y Breaks' if lb == refined_line_breaks[0] else "")
# Add the rightstrokes and leftstrokes based on X coordinate analysis to the plot
plt.scatter(rightstrokes, stagepos_y[rightstrokes], color='red', marker='x', label='Rightstrokes (X Analysis)')
plt.scatter(leftstrokes, stagepos_y[leftstrokes], color='blue', marker='o', label='Leftstrokes (X Analysis)')
plt.title(f'Stage Position Y with Line Breaks (Median Step Size: {median_step_size:.2f})')
plt.xlabel('Frame Number')
plt.ylabel('Stage Position Y')
plt.legend()
plt.show()
[docs]
def find_closest_stroke(frame, strokes):
# Find the closest frame in strokes to the given frame
closest_stroke = min(strokes, key=lambda x: abs(x - frame))
return closest_stroke
[docs]
def write_combined_stroke_info(input_path, rightstrokes, leftstrokes, refined_line_breaks):
with h5py.File(input_path, 'r+') as file:
# Get the length of the dataset
data_length = file['entry/data/stagepos_x_refined'].shape[0]
# Initialize arrays for linebreaks and streak directions
combined_lineinfo = np.zeros(data_length, dtype=int)
for index in range(data_length):
if index in refined_line_breaks:
closest_right = find_closest_stroke(index, rightstrokes)
closest_left = find_closest_stroke(index, leftstrokes)
if abs(closest_right - index) < abs(closest_left - index):
combined_lineinfo[index] = 1 # Closer to a right stroke
else:
combined_lineinfo[index] = -1 # Closer to a left stroke
else:
if index > 0:
combined_lineinfo[index] = combined_lineinfo[index - 1] # Carry forward the previous value
for i in combined_lineinfo:
log_print(combined_lineinfo[i])
# Write arrays to the HDF file
if 'entry/data/streakdirection' in file:
del file['entry/data/streakdirection']
file.create_dataset('entry/data/streakdirection', data=combined_lineinfo)
[docs]
def find_line_breaks(input_path):
rightstrokes, leftstrokes = find_strokes_by_xpos(input_path)
line_breaks, refined_line_breaks, median_step_size = find_strokes_by_ypos(input_path)
unmatched_leftstrokes, unmatched_rightstrokes = compare_strokes(leftstrokes, rightstrokes, refined_line_breaks)
if len(unmatched_leftstrokes) == 0 and len(unmatched_rightstrokes) == 0:
write_combined_stroke_info(input_path, rightstrokes, leftstrokes, refined_line_breaks)
else:
plot_stagepos_y_with_breaks(input_path)
raise ValueError("Unmatched strokes detected, check comaprison and ignore windows")
[docs]
def plot_stagepos_y_with_breaks_debug(input_path):
# Plotting line breaks based on y position for debugging
rightstrokes, leftstrokes = find_strokes_by_xpos(input_path)
line_breaks, refined_line_breaks, median_step_size = find_strokes_by_ypos(input_path)
with h5py.File(input_path, 'r') as file:
stagepos_y = file['entry/data/stagepos_y'][:]
streakdirection = file['entry/data/streakdirection'][:]
log_print(f'y-method: {len(refined_line_breaks)}')
plt.figure(figsize=(14, 7))
plt.plot(stagepos_y, label='Stage Position Y')
# Plot the initial line breaks based on Y coordinate analysis
#for lb in line_breaks:
# plt.axvline(x=lb, color='orange', linestyle=':', label='Initial Y Breaks' if lb == line_breaks[0] else "")
unmatched_leftstrokes, unmatched_rightstrokes = compare_strokes(leftstrokes, rightstrokes, refined_line_breaks)
log_print("Unmatched leftstrokes (X method):", unmatched_leftstrokes)
log_print("Unmatched rightstrokes (X method):", unmatched_rightstrokes)
# Plot the refined line breaks based on Y coordinate analysis
for lb in refined_line_breaks:
plt.axvline(x=lb, color='green', linestyle='--', label='Refined Y Breaks' if lb == refined_line_breaks[0] else "")
log_print("1")
# Add the rightstrokes and leftstrokes based on X coordinate analysis to the plot
plt.scatter(rightstrokes, stagepos_y[rightstrokes], color='red', marker='x', label='Rightstrokes (X Analysis)')
plt.scatter(leftstrokes, stagepos_y[leftstrokes], color='blue', marker='o', label='Leftstrokes (X Analysis)')
log_print("2")
# Scatter plot of stagepos_y values with colors based on streakdirection
for i in range(len(stagepos_y)):
color = 'red' if streakdirection[i] == 1 else 'purple' if streakdirection[i] == -1 else 'gray'
plt.scatter(i, stagepos_y[i], color=color)
log_print("3")
plt.title(f'Stage Position Y with Line Breaks (Median Step Size: {median_step_size:.2f})')
plt.xlabel('Frame Number')
plt.ylabel('Stage Position Y')
plt.legend()
# Save the plot as a PDF
plt.savefig(r"D:\20231204_Lys_Scilife\lys_2M_minus45deg_run_2023-12-04_16-50-53\plot", format='pdf')
plt.show()
[docs]
def get_three_line_bounds(input_path, frame_index):
with h5py.File(input_path, 'r') as workingfile:
streakdirection = workingfile['entry/data/streakdirection'][:]
# Identifying the bounds of the middle line
middle_line_bounds = get_line_bounds_from_index(streakdirection, frame_index)
# Identifying the bounds of the previous and next lines
prev_line_bounds = get_line_bounds_from_index(streakdirection, middle_line_bounds[0] - 1) if middle_line_bounds[0] > 0 else (None, None)
next_line_bounds = get_line_bounds_from_index(streakdirection, middle_line_bounds[1] + 1) if middle_line_bounds[1] < len(streakdirection) - 1 else (None, None)
return prev_line_bounds, middle_line_bounds, next_line_bounds
[docs]
def get_line_bounds_from_index(streakdirection, index):
current_direction = streakdirection[index]
start_index = index
while start_index > 0 and streakdirection[start_index - 1] == current_direction:
start_index -= 1
end_index = index
while end_index < len(streakdirection) - 1 and streakdirection[end_index + 1] == current_direction:
end_index += 1
return (start_index, end_index)
[docs]
def cross_correlate_lines(input_path, line_bounds1, line_bounds2, data_key):
with h5py.File(input_path, 'r') as workingfile:
data = workingfile[data_key][:]
line1 = data[line_bounds1[0]:line_bounds1[1]+1]
line2 = data[line_bounds2[0]:line_bounds2[1]+1]
correlation = np.correlate(line1, line2, mode='full')
return correlation
[docs]
def plot_correlation(correlation, title):
plt.figure(figsize=(10, 4))
plt.plot(correlation)
plt.title(title)
plt.xlabel('Lag')
plt.ylabel('Cross-correlation')
plt.grid(True)
plt.show()
# Example usage
#input_path = r'D:\20231204_Lys_Scilife\lys_2M_minus45deg_run_2023-12-04_16-50-53\lys_2M_minus45deg_run_2023-12-04_16-50-53.h5'
"""
data_key = 'entry/data/mean_intensities' # Replace with the key for your data
frame_index = 10000 # Replace with an index somewhere in your middle line
# Get bounds of the three lines
prev_line_bounds, middle_line_bounds, next_line_bounds = get_three_line_bounds(input_path, frame_index)
# Perform cross-correlation
correlation_prev_middle = cross_correlate_lines(input_path, prev_line_bounds, middle_line_bounds, data_key)
correlation_middle_next = cross_correlate_lines(input_path, middle_line_bounds, next_line_bounds, data_key)
# Plot the results
plot_correlation(correlation_prev_middle, "Cross-correlation: Previous Line vs Middle Line")
plot_correlation(correlation_middle_next, "Cross-correlation: Middle Line vs Next Line")
"""
[docs]
def moving_average(data, window_size):
""" Smooth data by performing a moving average. """
return np.convolve(data, np.ones(window_size)/window_size, mode='valid')
[docs]
def find_closest_minimum(x_data, y_data, target_x):
""" Find the closest minimum in y_data to the target_x in x_data """
# Find local minima
minima_indices, _ = find_peaks(-y_data)
if len(minima_indices) == 0:
return None, None # No minima found
# Find the closest minimum to target_x
closest_index = minima_indices[np.argmin(np.abs(x_data[minima_indices] - target_x))]
return x_data[closest_index], y_data[closest_index]
[docs]
def plot_line_profiles(input_path, index, window_size=5):
with h5py.File(input_path, 'r') as workingfile:
streakdirection = workingfile['entry/data/streakdirection'][:]
mean_intensities = workingfile['entry/data/mean_intensities'][:]
stagepos_x_refined = workingfile['entry/data/stagepos_x_refined'][:]
# Get bounds of the current, previous, and next lines
current_line_bounds = get_line_bounds_from_index(streakdirection, index)
prev_line_bounds = get_line_bounds_from_index(streakdirection, current_line_bounds[0] - 1) if current_line_bounds[0] > 0 else (None, None)
next_line_bounds = get_line_bounds_from_index(streakdirection, current_line_bounds[1] + 1) if current_line_bounds[1] < len(streakdirection) - 1 else (None, None)
plt.figure(figsize=(12, 6))
# Plot and process each line
for line_bounds, line_label in [(prev_line_bounds, 'Previous Line'),
(current_line_bounds, 'Current Line'),
(next_line_bounds, 'Next Line')]:
if line_bounds != (None, None):
x_line = stagepos_x_refined[line_bounds[0]:line_bounds[1] + 1]
y_line = mean_intensities[line_bounds[0]:line_bounds[1] + 1]
y_line_smooth = moving_average(y_line, window_size)
plt.plot(x_line[window_size - 1:], y_line_smooth, label=line_label) # Adjust x-axis to match smoothed data
# Find and annotate minima in the previous line
if prev_line_bounds != (None, None):
x_prev = stagepos_x_refined[prev_line_bounds[0]:prev_line_bounds[1] + 1]
y_prev = mean_intensities[prev_line_bounds[0]:prev_line_bounds[1] + 1]
y_prev_smooth = moving_average(y_prev, window_size)
min_index_prev = np.argmin(y_prev_smooth)
min_value_prev = y_prev_smooth[min_index_prev]
min_x_coordinate_prev = x_prev[min_index_prev + window_size - 1] # Adjust for window size
plt.scatter(min_x_coordinate_prev, min_value_prev, color='blue')
plt.annotate(f'Prev Min: {min_value_prev:.2f}\nX: {min_x_coordinate_prev:.2f}',
(min_x_coordinate_prev, min_value_prev),
textcoords="offset points",
xytext=(10, -10),
ha='center',
color='blue')
# Find and annotate closest minima in current and next lines
for line_bounds, line_color in [(current_line_bounds, 'orange'), (next_line_bounds, 'green')]:
if line_bounds != (None, None):
x_line = stagepos_x_refined[line_bounds[0]:line_bounds[1] + 1]
y_line = mean_intensities[line_bounds[0]:line_bounds[1] + 1]
y_line_smooth = moving_average(y_line, window_size)
x_line_smooth = x_line[window_size - 1:] # Adjust x-axis to match smoothed data
closest_min_x, closest_min_y = find_closest_minimum(x_line_smooth, y_line_smooth, min_x_coordinate_prev)
if closest_min_x is not None:
plt.scatter(closest_min_x, closest_min_y, color=line_color)
plt.annotate(f'Closest Min: {closest_min_y:.2f}\nX: {closest_min_x:.2f}',
(closest_min_x, closest_min_y),
textcoords="offset points",
xytext=(10, -10),
ha='center',
color=line_color)
plt.xlabel('Stage Position X (Refined)')
plt.ylabel('Mean Intensities')
plt.title('Line Profiles (Smoothed)')
plt.legend()
plt.show()
#index = 15000 # Replace with an index in your desired line
#plot_line_profiles(input_path, index)