from coseda.logging_utils import log_print
import h5py
from PIL import Image, ImageDraw, ImageEnhance
from sklearn.cluster import DBSCAN
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.colors as mcolors
import h5py
import cv2
import numpy as np
import matplotlib.cm as cm
from scipy.ndimage import gaussian_filter
import math
from collections import defaultdict
[docs]
def plot_map_peaks(h5file_path, zdim, on_pick_event_handler=None):
with h5py.File(h5file_path, 'r') as workingfile:
zpath = 'entry/data/' + str(zdim)
if 'entry/data/stagepos_x_refined' in workingfile:
stagepos_x = workingfile['entry/data/stagepos_x_refined'][:]
else:
stagepos_x = workingfile['entry/data/stagepos_x'][:]
if 'entry/data/stagepos_y_refined' in workingfile:
stagepos_y = workingfile['entry/data/stagepos_y_refined'][:]
else:
stagepos_x = workingfile['entry/data/stagepos_y'][:]
zdimension = workingfile[zpath][:]
valid_mask = (stagepos_x != 999999) & (stagepos_y != 999999)
# Apply the mask to stage positions and zdimension
stagepos_x_filtered = stagepos_x[valid_mask]
stagepos_y_filtered = stagepos_y[valid_mask]
zdimension_filtered = zdimension[valid_mask]
# Determine the scaling and translation factors for the coordinates
x_min, x_max = stagepos_x_filtered.min(), stagepos_x_filtered.max()
y_min, y_max = stagepos_y_filtered.min(), stagepos_y_filtered.max()
log_print('x range: ' + str(x_max - x_min))
log_print(x_min, x_max)
log_print('y range: ' + str(y_max - y_min))
log_print(y_min, y_max)
# Normalize and visualize
plt.figure()
plt.scatter(stagepos_x_filtered, stagepos_y_filtered, c=zdimension_filtered, cmap='inferno_r', marker='.', s=100)
#plt.colorbar(label='Number of Peaks')
plt.title('Reconstructed Image')
plt.axis('equal') # Ensure equal scaling of the x and y axes
plt.savefig('intensity_map.png')
#plt.show()
scatter = plt.scatter(stagepos_x_filtered, stagepos_y_filtered, c=zdimension_filtered, cmap='inferno_r', marker='.', s=200, picker=True)
plt.colorbar(scatter, label='Number of Peaks')
if on_pick_event_handler is not None:
# Connect the pick event handler
plt.gcf().canvas.mpl_connect('pick_event', on_pick_event_handler)
plt.show()
[docs]
def plot_crystamorphusvac(input_path, offset = None):
with h5py.File(input_path, 'r') as workingfile:
stagepos_x = workingfile['entry/data/stagepos_x_refined'][:]
stagepos_y = workingfile['entry/data/stagepos_y'][:]
total_intensities = workingfile['entry/data/mean_intensities'][:]
npeaks = workingfile['entry/data/nPeaks'][:]
streakdirection = workingfile['entry/data/streakdirection'][:]
if offset is not None:
# Apply offset to stagepos_x where streakdirection is -1
stagepos_x[streakdirection == -1] += offset
# Determine the scaling and translation factors for the coordinates
x_min, x_max = stagepos_x.min(), stagepos_x.max()
y_min, y_max = stagepos_y.min(), stagepos_y.max()
# print('x range: ' + str(x_max - x_min))
# print(x_min, x_max)
# print('y range: ' + str(y_max - y_min))
# print(y_min, y_max)
flag = [None] * len(total_intensities)
for i in range(len(total_intensities)):
if total_intensities[i] >= 0.1 * total_intensities.max() and npeaks[i] > 15:
flag[i] = 'cryst' # Mark as crystalline
elif total_intensities[i] >= 0.5 * total_intensities.max():
flag[i] = 'amorph' # Mark as amorphous
elif total_intensities[i] >= 0.10 * total_intensities.max():
flag[i] = 'vac' # Mark as amorphous
else:
flag[i] = 'copper' # Mark as amorphous
# If neither condition is met, flag[i] remains None
color_map = {'cryst': 'magenta', 'amorph': 'blue', 'vac': 'lightsteelblue', 'copper': 'sandybrown'}
colors = [color_map[flag_value] for flag_value in flag]
# Normalize and visualize
plt.figure()
scatter = plt.scatter(stagepos_x, stagepos_y, c=colors, marker='.', s=1)
# Create a legend
legend_labels = {'cryst': 'Crystalline', 'amorph': 'Amorphous', 'vac': 'Vacuum/Copper', None: 'Copper'}
patches = [mpatches.Patch(color=color, label=label) for label, color in color_map.items()]
plt.legend(handles=patches, title="Classification")
plt.title('Flag Mapping')
plt.axis('equal') # Ensure equal scaling of the x and y axes
plt.savefig('flag_map.png')
plt.show()
[docs]
def plot_crystamorphusvac_overlap(input_path):
with h5py.File(input_path, 'r') as file:
stagepos_x = file['entry/data/stagepos_x_refined'][:]
stagepos_y = file['entry/data/stagepos_y'][:]
total_intensities = file['entry/data/mean_intensities'][:]
npeaks = file['entry/data/nPeaks'][:]
# Determine the scaling and translation factors for the coordinates
x_min, x_max = stagepos_x.min(), stagepos_x.max()
y_min, y_max = stagepos_y.min(), stagepos_y.max()
# Generate flag array outside the file context
flag = [None] * len(total_intensities)
for i in range(len(total_intensities)):
if total_intensities[i] >= 0.1 * total_intensities.max() and npeaks[i] > 20:
flag[i] = 'cryst' # Mark as crystalline
elif total_intensities[i] >= 0.3 * total_intensities.max():
flag[i] = 'amorph' # Mark as amorphous
#elif total_intensities[i] >= 0.05 * total_intensities.max():
# flag[i] = 'vac/copper' # Mark as amorphous
else:
flag[i] = 'vac' # Mark as amorphous
# If neither condition is met, flag[i] remains None
# Interpolation code...
overlap_extent = 0.5 # Assuming 50% overlap
x_interpolated = np.linspace(x_min, x_max, int(len(stagepos_x) / overlap_extent))
y_interpolated = np.linspace(y_min, y_max, int(len(stagepos_y) / overlap_extent))
interpolated_flag = np.empty((len(x_interpolated), len(y_interpolated)), dtype=object)
for i, x_val in enumerate(x_interpolated):
for j, y_val in enumerate(y_interpolated):
# Find the closest original pixel
dist = np.sqrt((stagepos_x - x_val)**2 + (stagepos_y - y_val)**2)
closest_idx = np.argmin(dist)
interpolated_flag[i, j] = flag[closest_idx]
color_map = {'cryst': 'magenta', 'amorph': 'blue', 'vac': 'lightsteelblue', None: 'sandybrown'}
# Create a color array for the interpolated grid
interpolated_colors = np.vectorize(color_map.get)(interpolated_flag)
# Normalize and visualize
plt.figure()
plt.imshow(interpolated_colors, origin='lower', extent=[x_min, x_max, y_min, y_max], aspect='auto')
# Create a legend
legend_labels = {'cryst': 'Crystalline', 'amorph': 'Amorphous', 'vac': 'Vacuum/Copper', None: 'Copper'}
patches = [mpatches.Patch(color=color, label=label) for label, color in color_map.items()]
plt.legend(handles=patches, title="Classification")
plt.title('Interpolated Flag Mapping')
plt.axis('equal') # Ensure equal scaling of the x and y axes
plt.savefig('interpolated_flag_map.png')
plt.show()
"""
def linearize_map(input_path, threshold, beam_diameter):
with h5py.File(input_path, 'r') as file:
stagepos_x = file['entry/data/stagepos_x_refined'][:]
stagepos_y = file['entry/data/stagepos_y'][:]
nPeaks = file['entry/data/nPeaks'][:]
segments = []
start_index = None
y_sum = 0
# Identify segments where nPeaks exceeds the threshold
for i in range(len(nPeaks)):
if nPeaks[i] > threshold:
if start_index is None:
start_index = i
y_sum += stagepos_y[i]
else:
if start_index is not None:
segment_length = i - start_index
avg_y = y_sum / segment_length if segment_length > 0 else 0
startx = stagepos_x[start_index] - beam_diameter / 2
endx = stagepos_x[i - 1] + beam_diameter / 2
segments.append(
(startx, stagepos_y[start_index], endx, stagepos_y[i - 1], avg_y, start_index, i - 1)
)
start_index = None
y_sum = 0
# Check for the case where the last frame is above the threshold
if start_index is not None:
segment_length = len(nPeaks) - start_index
avg_y = y_sum / segment_length if segment_length > 0 else 0
startx = stagepos_x[start_index] - beam_diameter / 2
endx = stagepos_x[-1] + beam_diameter / 2
segments.append(
(startx, stagepos_y[start_index], endx, stagepos_y[-1], avg_y, start_index, len(nPeaks) - 1)
)
# Merge overlapping segments within a line
merged_segments = []
i = 0
while i < len(segments):
current_seg = segments[i]
j = i + 1
while j < len(segments):
next_seg = segments[j]
if current_seg[6] == next_seg[6] and current_seg[2] >= next_seg[0]:
merged_startx = min(current_seg[0], next_seg[0])
merged_endx = max(current_seg[2], next_seg[2])
merged_start_index = min(current_seg[5], next_seg[5])
merged_end_index = max(current_seg[6], next_seg[6])
current_seg = (
merged_startx,
current_seg[1],
merged_endx,
current_seg[3],
current_seg[4],
merged_start_index,
merged_end_index,
)
j += 1
else:
break
merged_segments.append(current_seg)
i = j
segments = merged_segments
# Plotting the intensity map
plt.figure()
plt.scatter(stagepos_x, stagepos_y, c='blue', marker='.', s=200)
for segment in segments:
plt.plot([segment[0], segment[2]], [segment[1], segment[3]], 'r-')
plt.title(f'Intensity Mapping with Average Y Segments (Frames >{threshold} Peaks)')
plt.axis('equal')
plt.show()
# Calculate y distances between consecutive segments
y_distances = [abs(segments[i + 1][1] - segments[i][3]) for i in range(len(segments) - 1)]
# Plotting the histogram of y distances
plt.figure()
plt.hist(y_distances, bins=30, color='green', edgecolor='black')
plt.title('Histogram of Y Distances Between Consecutive Segments')
plt.xlabel('Y Distance')
plt.ylabel('Frequency')
plt.show()
# Use local standard deviation to set a dynamic threshold for each segment
local_std_threshold = [
np.std(y_distances[max(0, i - 5):min(len(y_distances), i + 5)])
for i in range(len(y_distances))
]
mean_local_std = np.mean(local_std_threshold)
std_local_std = np.std(local_std_threshold)
dynamic_threshold = [mean_local_std + 2 * std_local_std for _ in range(len(y_distances))]
# Classify as line break or inaccuracy
line_breaks = [y_distances[i] > dynamic_threshold[i] for i in range(len(y_distances))]
line_number = 1
segment_index = 1
for i in range(len(segments)):
if i < len(line_breaks) and line_breaks[i]:
line_number += 1
segments[i] = segments[i] + (line_number, segment_index)
segment_index += 1
# Plotting the intensity map with line numbers
plt.figure()
plt.scatter(stagepos_x, stagepos_y, c='blue', marker='.', s=200)
for segment in segments:
plt.plot([segment[0], segment[2]], [segment[7], segment[7]], 'r-')
plt.text((segment[0] + segment[2]) / 2, segment[7], str(segment[8]))
plt.title(f'Mapping segments based on assigned line numbers (Frames >{threshold} Peaks)')
plt.xlabel('Stagepos X')
plt.ylabel('Line Number')
plt.gca().invert_yaxis()
plt.show()
return segments
def is_point_near_line(point, line_start, line_end, tolerance):
# Function to check if a point is near a given line within a certain tolerance
# Using the distance formula from point to line
x0, y0 = point
x1, y1 = line_start
x2, y2 = line_end
dist = abs((y2 - y1)*x0 - (x2 - x1)*y0 + x2*y1 - y2*x1) / np.sqrt((y2 - y1)**2 + (x2 - x1)**2)
return dist <= tolerance
def find_best_matching_lines(points, tolerance):
best_lines = []
max_count = 0
for i in range(len(points)):
for j in range(i+1, len(points)):
current_count = 0
line_start = points[i]
line_end = points[j]
for point in points:
if is_point_near_line(point, line_start, line_end, tolerance):
current_count += 1
if current_count > max_count:
best_lines = [(line_start, line_end)]
max_count = current_count
elif current_count == max_count:
best_lines.append((line_start, line_end))
return best_lines, max_count
def find_best_matching_lines_without_duplicates(points, tolerance, num_lines=3):
lines_with_points = []
for i in range(len(points)):
for j in range(i+1, len(points)):
line_start = points[i]
line_end = points[j]
line_points = []
for point in points:
if is_point_near_line(point, line_start, line_end, tolerance):
line_points.append(point)
# Add the line only if it has 3 or more points close to it
if len(line_points) >= 4:
lines_with_points.append((line_start, line_end, line_points))
# Sort lines based on the number of points near them
lines_with_points.sort(key=lambda x: len(x[2]), reverse=True)
best_lines = []
used_points = set()
for line in lines_with_points:
if len(best_lines) >= num_lines:
break
line_start, line_end, line_points = line
if not any(p in used_points for p in line_points):
best_lines.append((line_start, line_end))
used_points.update(line_points)
return best_lines
def get_extended_line_points(line_start, line_end, img_shape):
Extend the line to the image boundaries and return the new line start and end points.
# Unpack image shape and line points
img_height, img_width = img_shape
x0, y0 = line_start
x1, y1 = line_end
# Calculate line coefficients: y = mx + b
if x1 - x0 != 0:
m = (y1 - y0) / (x1 - x0)
b = y0 - m * x0
# Find intersection with image boundaries
y_start = 0
x_start = -b / m if m != 0 else x0
y_end = img_height
x_end = (y_end - b) / m if m != 0 else x0
else:
# Vertical line
x_start, x_end = x0, x0
y_start, y_end = 0, img_height
# Ensure points are within image boundaries
x_start, x_end = max(0, min(x_start, img_width)), max(0, min(x_end, img_width))
y_start, y_end = max(0, min(y_start, img_height)), max(0, min(y_end, img_height))
return (x_start, y_start), (x_end, y_end)
def calculate_angle(line_start, line_end):
Calculate the angle (in degrees) of the line with respect to the x-axis.
x0, y0 = line_start
x1, y1 = line_end
angle_rad = np.arctan2(y1 - y0, x1 - x0)
angle_deg = np.degrees(angle_rad)
if angle_deg > 180:
angle_deg = 360 - angle_deg
elif angle_deg < 0:
angle_deg = -angle_deg
return angle_deg
def identify_crystals(input_path, frame_indices, gamma, tolerance, num_lines):
with h5py.File(input_path, 'r') as file:
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for i, frame_index in enumerate(frame_indices):
image_data = file['entry/data/images'][frame_index, :, :]
image_shape = image_data.shape
image_data_normalized = (image_data - image_data.min()) / (image_data.max() - image_data.min())
image_data_gamma_corrected = np.power(image_data_normalized, gamma)
colored_image = cm.get_cmap('inferno')(image_data_gamma_corrected)
axes[i].imshow(colored_image, aspect='auto')
peak_x_positions = file['entry/data/peakXPosRaw'][frame_index, :]
peak_y_positions = file['entry/data/peakYPosRaw'][frame_index, :]
points = [(x, y) for x, y in zip(peak_x_positions, peak_y_positions) if (x, y) != (0, 0)]
best_lines = find_best_matching_lines_without_duplicates(points, tolerance, num_lines)
axes[i].scatter(peak_x_positions, peak_y_positions, s=50, facecolors='none', edgecolors='white')
for line_start, line_end in best_lines:
extended_start, extended_end = get_extended_line_points(line_start, line_end, image_shape)
angle = calculate_angle(line_start, line_end)
axes[i].plot([extended_start[0], extended_end[0]], [extended_start[1], extended_end[1]], 'r-')
axes[i].text(line_start[0], line_start[1], f'{angle:.1f}°', color='red', fontsize=10)
axes[i].set_title(f'Frame {frame_index}')
axes[i].axis('equal')
plt.show()
# frame = 22915
# frame_indices = [frame-1, frame, frame+1]
# gamma_value = 0.4
# tolerance = 8
# num_lines = 10
# identify_crystals(input_path,frame_indices,gamma_value,tolerance, num_lines)
def find_clusters(input_path, distance_threshold=0.00000005):
with h5py.File(input_path, 'r') as file:
stagepos_x = file['entry/data/stagepos_x_refined'][:]
stagepos_y = file['entry/data/stagepos_y'][:]
nPeaks = file['entry/data/nPeaks'][:]
# Filter frames with nPeaks > 10
filtered_indices = np.where(nPeaks > 10)[0]
filtered_x = stagepos_x[filtered_indices]
filtered_y = stagepos_y[filtered_indices]
# Initialize clusters
clusters = [-1] * len(filtered_indices) # -1 indicates no cluster
current_cluster = 0
for i in range(len(filtered_indices)):
if clusters[i] == -1: # not yet assigned to a cluster
clusters[i] = current_cluster
for j in range(i + 1, len(filtered_indices)):
if clusters[j] == -1: # also not yet assigned
distance = np.sqrt((filtered_x[i] - filtered_x[j])**2 + (filtered_y[i] - filtered_y[j])**2)
if distance <= distance_threshold:
clusters[j] = current_cluster
current_cluster += 1
# Visualize the clusters
plt.figure()
for i, index in enumerate(filtered_indices):
plt.scatter(stagepos_x[index], stagepos_y[index], label=f'Cluster {clusters[i]}')
plt.title('Intensity Mapping with Clusters')
plt.axis('equal')
plt.legend()
plt.show()
def find_crystals(input_path):
with h5py.File(input_path, 'r') as file:
stagepos_x = file['entry/data/stagepos_x_refined'][:]
stagepos_y = file['entry/data/stagepos_y'][:]
nPeaks = file['entry/data/nPeaks'][:]
# Determine the scaling and translation factors for the coordinates
x_min, x_max = stagepos_x.min(), stagepos_x.max()
y_min, y_max = stagepos_y.min(), stagepos_y.max()
log_print('x range: ' + str(x_max - x_min))
log_print(x_min, x_max)
log_print('y range: ' + str(y_max - y_min))
log_print(y_min, y_max)
# Normalize and visualize
plt.figure()
scatter = plt.scatter(stagepos_x, stagepos_y, c=nPeaks, cmap='inferno_r', marker='.', s=200)
plt.colorbar(label='Normalized Total Intensity')
plt.title('Intensity Mapping')
plt.axis('equal')
# Add hover effect
cursor = mplcursors.cursor(scatter, hover=True)
cursor.connect("add", lambda sel: sel.annotation.set_text(
f'index: {sel.target.index}, nPeaks: {nPeaks[sel.target.index]:.2f}'))
plt.savefig('intensity_map.png')
plt.show()
def find_crystals(input_path, threshold):
with h5py.File(input_path, 'r') as file:
stagepos_x = file['entry/data/stagepos_x_refined'][:]
stagepos_y = file['entry/data/stagepos_y'][:]
nPeaks = file['entry/data/nPeaks'][:]
# Determine the scaling and translation factors for the coordinates
x_min, x_max = stagepos_x.min(), stagepos_x.max()
y_min, y_max = stagepos_y.min(), stagepos_y.max()
# Normalize and visualize
plt.figure()
colors = ['red' if value > threshold else 'blue' for value in nPeaks] # Highlighting condition
scatter = plt.scatter(stagepos_x, stagepos_y, c=colors, marker='.', s=200)
plt.colorbar(label='Normalized Total Intensity')
plt.title('Intensity Mapping')
plt.axis('equal')
# Add hover effect
cursor = mplcursors.cursor(scatter, hover=True)
cursor.connect("add", lambda sel: sel.annotation.set_text(
f'index: {sel.target.index}, nPeaks: {nPeaks[sel.target.index]:.2f}'))
plt.savefig('intensity_map.png')
plt.show()
def linearize_map(input_path, threshold, beam_diameter):
with h5py.File(input_path, 'r') as file:
stagepos_x = file['entry/data/stagepos_x_refined'][:]
stagepos_y = file['entry/data/stagepos_y'][:]
nPeaks = file['entry/data/nPeaks'][:]
segments = []
start_index = None
y_sum = 0
# Identify segments where nPeaks exceeds the threshold
for i in range(len(nPeaks)):
if nPeaks[i] > threshold:
if start_index is None:
start_index = i
y_sum += stagepos_y[i]
else:
if start_index is not None:
segment_length = i - start_index
avg_y = y_sum / segment_length if segment_length > 0 else 0
startx = stagepos_x[start_index] - beam_diameter / 2
endx = stagepos_x[i-1] + beam_diameter / 2
segments.append((startx, stagepos_y[start_index], endx, stagepos_y[i-1], avg_y, start_index, i-1))
start_index = None
y_sum = 0
# Check for the case where the last frame is above the threshold
if start_index is not None:
segment_length = len(nPeaks) - start_index
avg_y = y_sum / segment_length if segment_length > 0 else 0
startx = stagepos_x[start_index] - beam_diameter / 2
endx = stagepos_x[-1] + beam_diameter / 2
segments.append((startx, stagepos_y[start_index], endx, stagepos_y[-1], avg_y, start_index, len(nPeaks) - 1))
# Merge overlapping segments within a line:
merged_segments = []
i = 0
while i < len(segments):
# Current segment to compare with others
current_seg = segments[i]
# Iterate over the next segments to check for overlap
j = i + 1
while j < len(segments):
next_seg = segments[j]
# Check if the segments are in the same line and overlap
if current_seg[6] == next_seg[6] and current_seg[2] >= next_seg[0]:
# Merge segments: use smaller start value, larger end value, and update indices
merged_startx = min(current_seg[0], next_seg[0])
merged_endx = max(current_seg[2], next_seg[2])
merged_start_index = min(current_seg[5], next_seg[5])
merged_end_index = max(current_seg[6], next_seg[6])
# Create new merged segment
current_seg = (merged_startx, current_seg[1], merged_endx, current_seg[3], current_seg[4], merged_start_index, merged_end_index, current_seg[7], current_seg[8], current_seg[9])
# Skip the next_seg as it's now merged
j += 1
else:
break
# Add the current (potentially merged) segment to the result
merged_segments.append(current_seg)
i = j # Move to the next non-overlapping segment
segments = merged_segments
# Plotting the intensity map
plt.figure()
plt.scatter(stagepos_x, stagepos_y, c='blue', marker='.', s=200)
# Plot lines for each segment using average y-coordinate
for segment in segments:
plt.plot([segment[0], segment[2]], [segment[1], segment[3]], 'r-')
#plt.plot([segment[0], segment[2]], [segment[4], segment[4]], 'g-')
plt.title(f'Intensity Mapping with Average Y Segments (Frames >{threshold} Peaks)')
plt.axis('equal')
plt.show()
# Calculate y distances between consecutive segments
y_distances = [abs(segments[i+1][1] - segments[i][3]) for i in range(len(segments) - 1)]
# Plotting the histogram of y distances
plt.figure()
plt.hist(y_distances, bins=30, color='green', edgecolor='black') # Adjust the number of bins as needed
plt.title('Histogram of Y Distances Between Consecutive Segments')
plt.xlabel('Y Distance')
plt.ylabel('Frequency')
plt.show()
# Use local standard deviation to set a dynamic threshold for each segment
local_std_threshold = [np.std(y_distances[max(0, i-5):min(len(y_distances), i+5)]) for i in range(len(y_distances))]
mean_local_std = np.mean(local_std_threshold)
std_local_std = np.std(local_std_threshold)
dynamic_threshold = [mean_local_std + 2 * std_local_std for _ in range(len(y_distances))]
# Classify as line break or inaccuracy
line_breaks = [y_distances[i] > dynamic_threshold[i] for i in range(len(y_distances))]
line_number = 1
segment_index = 1
for i in range(len(segments)):
# Update line number when a line break is detected
if i < len(line_breaks) and line_breaks[i]:
line_number += 1
# Create a new tuple with the existing segment data and the line number
segments[i] = segments[i] + (line_number, segment_index)
segment_index += 1
# Plotting the intensity map with lines numbers
plt.figure()
plt.scatter(stagepos_x, stagepos_y, c='blue', marker='.', s=200)
# Plot lines for each segment using average y-coordinate
for segment in segments:
plt.plot([segment[0], segment[2]], [segment[7], segment[7]], 'r-')
plt.text((segment[0] + segment[2]) / 2, segment[7], str(segment[8]))
plt.title(f'Mapping segments based on assigned line numbers (Frames >{threshold} Peaks)')
plt.xlabel('Stagepos X')
plt.ylabel('Line Number')
plt.gca().invert_yaxis()
#plt.axis('equal')
plt.show()
'''
# Plotting the histogram with dynamic thresholds
plt.figure()
colors = ['red' if line_breaks[i] else 'blue' for i in range(len(y_distances))]
plt.bar(range(len(y_distances)), y_distances, color=colors)
plt.title('Y Distance Between Consecutive Segments with Dynamic Line Breaks Highlighted')
plt.xlabel('Segment Index')
plt.ylabel('Y Distance')
plt.show()
'''
return segments
"""
[docs]
def linearize_map(input_path, threshold, beam_diameter):
with h5py.File(input_path, 'r') as file:
stagepos_x = file['entry/data/stagepos_x_refined'][:]
stagepos_y = file['entry/data/stagepos_y'][:]
nPeaks = file['entry/data/nPeaks'][:]
segments = []
start_index = None
y_sum = 0
# Identify segments where nPeaks exceeds the threshold
for i in range(len(nPeaks)):
if nPeaks[i] > threshold:
if start_index is None:
start_index = i
y_sum += stagepos_y[i]
else:
if start_index is not None:
segment_length = i - start_index
avg_y = y_sum / segment_length if segment_length > 0 else 0
startx = stagepos_x[start_index] - beam_diameter / 2
endx = stagepos_x[i - 1] + beam_diameter / 2
segments.append(
(startx, stagepos_y[start_index], endx, stagepos_y[i - 1], avg_y, start_index, i - 1)
)
start_index = None
y_sum = 0
# Check for the case where the last frame is above the threshold
if start_index is not None:
segment_length = len(nPeaks) - start_index
avg_y = y_sum / segment_length if segment_length > 0 else 0
startx = stagepos_x[start_index] - beam_diameter / 2
endx = stagepos_x[-1] + beam_diameter / 2
segments.append(
(startx, stagepos_y[start_index], endx, stagepos_y[-1], avg_y, start_index, len(nPeaks) - 1)
)
# Merge overlapping segments within a line
merged_segments = []
i = 0
while i < len(segments):
current_seg = segments[i]
j = i + 1
while j < len(segments):
next_seg = segments[j]
if current_seg[6] == next_seg[6] and current_seg[2] >= next_seg[0]:
merged_startx = min(current_seg[0], next_seg[0])
merged_endx = max(current_seg[2], next_seg[2])
merged_start_index = min(current_seg[5], next_seg[5])
merged_end_index = max(current_seg[6], next_seg[6])
current_seg = (
merged_startx,
current_seg[1],
merged_endx,
current_seg[3],
current_seg[4],
merged_start_index,
merged_end_index,
)
j += 1
else:
break
merged_segments.append(current_seg)
i = j
segments = merged_segments
# Plotting the intensity map
plt.figure()
plt.scatter(stagepos_x, stagepos_y, c='blue', marker='.', s=200)
for segment in segments:
plt.plot([segment[0], segment[2]], [segment[1], segment[3]], 'r-')
plt.title(f'Intensity Mapping with Average Y Segments (Frames >{threshold} Peaks)')
plt.axis('equal')
plt.show()
# Calculate y distances between consecutive segments
y_distances = [abs(segments[i + 1][1] - segments[i][3]) for i in range(len(segments) - 1)]
# Plotting the histogram of y distances
plt.figure()
plt.hist(y_distances, bins=30, color='green', edgecolor='black')
plt.title('Histogram of Y Distances Between Consecutive Segments')
plt.xlabel('Y Distance')
plt.ylabel('Frequency')
plt.show()
# Use local standard deviation to set a dynamic threshold for each segment
local_std_threshold = [
np.std(y_distances[max(0, i - 5):min(len(y_distances), i + 5)])
for i in range(len(y_distances))
]
mean_local_std = np.mean(local_std_threshold)
std_local_std = np.std(local_std_threshold)
dynamic_threshold = [mean_local_std + 2 * std_local_std for _ in range(len(y_distances))]
# Classify as line break or inaccuracy
line_breaks = [y_distances[i] > dynamic_threshold[i] for i in range(len(y_distances))]
line_number = 1
segment_index = 1
for i in range(len(segments)):
if i < len(line_breaks) and line_breaks[i]:
line_number += 1
segments[i] = segments[i] + (line_number, segment_index)
segment_index += 1
# Plotting the intensity map with line numbers
plt.figure()
plt.scatter(stagepos_x, stagepos_y, c='blue', marker='.', s=200)
for segment in segments:
plt.plot([segment[0], segment[2]], [segment[7], segment[7]], 'r-')
plt.text((segment[0] + segment[2]) / 2, segment[7], str(segment[8]))
plt.title(f'Mapping segments based on assigned line numbers (Frames >{threshold} Peaks)')
plt.xlabel('Stagepos X')
plt.ylabel('Line Number')
plt.gca().invert_yaxis()
plt.show()
return segments
[docs]
def cluster_segments(input_path, threshold, beam_diameter):
segments = linearize_map(input_path, threshold, beam_diameter)
# Group segments by line number
lines = {}
for segment in segments:
line_number = segment[7] # The line number is the 8th element in the segment tuple
lines.setdefault(line_number, []).append(segment)
# List to store tuples of overlapping segment indices
overlapping_segments = []
# # Iterate through each line and compare with the next line
# for line_number in sorted(lines.keys())[:-1]: # Exclude the last line
# current_line_segments = lines[line_number]
# next_line_segments = lines[line_number + 1]
# # Debugging: print the line being processed
# print(f"Processing line {line_number} and line {line_number + 1}")
# # Iterate over segments in the current line
# for seg_a in current_line_segments:
# # Iterate over segments in the next line
# for seg_b in next_line_segments:
# # Debugging: print segments being compared
# print(f"Comparing segments {seg_a[8]} and {seg_b[8]} with x-ranges {seg_a[0]}-{seg_a[2]} and {seg_b[0]}-{seg_b[2]}")
# # Check for x overlap
# if seg_a[2] >= seg_b[0] and seg_a[0] <= seg_b[2]:
# # If overlap, add tuple of segment indices to the list
# overlapping_segments.append((seg_a[8], seg_b[8])) # segment index
# Assuming 'segments' is a list of segments where each segment has the line number at index [7]
for i, seg_a in enumerate(segments):
line_a = seg_a[7] # Get the line number of seg_a
# Check against other segments
for j, seg_b in enumerate(segments):
if i != j: # Avoid comparing a segment with itself
line_b = seg_b[7]
# Check if seg_b is on the same line or adjacent line
if line_b in [line_a, line_a - 1, line_a + 1]:
# Compare x-ranges for overlap (adjust indices as per your data structure)
if seg_a[2] >= seg_b[0] and seg_a[0] <= seg_b[2]:
# Add overlapping segments to the list
overlapping_segments.append((i, j))
# overlapping_segments contains tuples of indices of overlapping segments between consecutive lines
#print(segments)
log_print(overlapping_segments)
log_print(f'Number of segments: {len(segments)}')
log_print(f'Number of overlapping segments: {len(overlapping_segments)}')
clusters = []
for overlapping_segment in overlapping_segments:
seg_a, seg_b = overlapping_segment
a_in_cluster = b_in_cluster = None
# Check if either of the segments already exists in any cluster
for i, cluster in enumerate(clusters):
if seg_a in cluster:
a_in_cluster = i
if seg_b in cluster:
b_in_cluster = i
if a_in_cluster is None and b_in_cluster is None:
# Neither segment is in any cluster, create a new cluster
clusters.append({seg_a, seg_b})
elif a_in_cluster is not None and b_in_cluster is None:
# Segment A is in a cluster, but B is not, add B to A's cluster
clusters[a_in_cluster].add(seg_b)
elif b_in_cluster is not None and a_in_cluster is None:
# Segment B is in a cluster, but A is not, add A to B's cluster
clusters[b_in_cluster].add(seg_a)
elif a_in_cluster != b_in_cluster:
# Both segments are in different clusters, merge these clusters
clusters[a_in_cluster].update(clusters[b_in_cluster])
del clusters[b_in_cluster]
def merge_overlapping_clusters(clusters):
i = 0
while i < len(clusters) - 1:
merged = False
for j in range(i + 1, len(clusters)):
if clusters[i].intersection(clusters[j]):
clusters[i].update(clusters[j])
del clusters[j]
merged = True
break
if not merged:
i += 1
return clusters
# Merge clusters with overlapping indices
clusters = merge_overlapping_clusters(clusters)
#print(clusters)
length1 = 0 # Cluster size 1
length2 = 0 # Cluster size 2
length3 = 0 # Cluster size 3
length4 = 0 # Cluster size 4
length5 = 0 # Cluster size 5
length6 = 0 # Cluster size 6
length7 = 0 # Cluster size 7
lengthlarger = 0 # Cluster size 8 or more
for cluster in clusters:
if len(cluster) == 1: # Cluster size 1
length1 += 1
elif len(cluster) == 2: # Cluster size 2
length2 += 1
elif len(cluster) == 3: # Cluster size 3
length3 += 1
elif len(cluster) == 4: # Cluster size 4
length4 += 1
elif len(cluster) == 5: # Cluster size 5
length5 += 1
elif len(cluster) == 6: # Cluster size 6
length6 += 1
elif len(cluster) == 7: # Cluster size 7
length7 += 1
else: # Cluster size 8 or more
lengthlarger += 1
# Print the count of each cluster size
log_print("Cluster size 1:", length1)
log_print("Cluster size 2:", length2)
log_print("Cluster size 3:", length3)
log_print("Cluster size 4:", length4)
log_print("Cluster size 5:", length5)
log_print("Cluster size 6:", length6)
log_print("Cluster size 7:", length7)
log_print("Cluster size 8+:", lengthlarger)
# Assuming 'segments' is your list of segments and 'clusters' is the list of sets of segment indices
segment_to_cluster = {}
# Create mapping from segment index to cluster index
for cluster_idx, cluster in enumerate(clusters):
for seg_index in cluster:
segment_to_cluster[seg_index] = cluster_idx
# Update segments with their respective cluster index or None
updated_segments = []
for segment in segments:
seg_index = segment[8] # Assuming the 9th element is the unique segment index
cluster_idx = segment_to_cluster.get(seg_index, None)
updated_segment = segment + (cluster_idx,) # Create a new tuple with the cluster index
updated_segments.append(updated_segment)
# Now, 'segments' contains the updated segments with cluster indices
#print(updated_segments)
# Print the clusters
#for cluster_index, cluster_segments in clusters.items():
# print(f"Cluster {cluster_index}: {cluster_segments}")
#for segment in updated_segments:
# print(f'{segment[8]} is part of {segment[9]}')
# Generate a color map with a unique color for each cluster
num_clusters = len(clusters)
colors = plt.cm.get_cmap('hsv', num_clusters)
for segment in updated_segments:
# Determine the cluster index of the segment
cluster_index = segment[9] # The new cluster index is now the 9th element
segment_index = segment[8]
line_number = segment[7]
color = colors(cluster_index) if cluster_index is not None else 'gray'
# Draw the segment line
plt.plot([segment[0], segment[2]], [segment[7], segment[7]], color=color, linewidth=4)
# Add text label for the cluster index at the center of the segment
mid_x = (segment[0] + segment[2]) / 2
plt.text(mid_x, segment[7] + 0.3, "seg"+ str(segment_index)+ "/cl" + str(cluster_index) + "/ln" + str(line_number), color=color, fontsize=8, ha='center', va='center')
plt.title(f'Mapping segments based on assigned line numbers (Frames >{threshold} Peaks)')
plt.xlabel('Stagepos X')
plt.ylabel('Line Number')
plt.gca().invert_yaxis()
plt.show()
[docs]
def plot_nPeaks_histogram(input_path):
with h5py.File(input_path, 'r') as file:
nPeaks = file['entry/data/nPeaks'][:]
# Plotting the histogram
plt.figure()
plt.hist(nPeaks, bins=5 , color='blue', edgecolor='black') # Adjust the number of bins as needed
plt.title('Histogram of nPeaks')
plt.xlabel('nPeaks Value')
plt.ylabel('Frequency')
plt.savefig('nPeaks_histogram.png')
plt.show()
[docs]
def find_matrix(input_path):
with h5py.File(input_path, 'r') as file:
stagepos_x = file['entry/data/stagepos_x_refined'][:]
stagepos_y = file['entry/data/stagepos_y'][:]
nPeaks = file['entry/data/nPeaks'][:]
xrange = max(stagepos_x) - min(stagepos_x)
turningzone = False
turncount = 0
for i in range(0, len(stagepos_x), 10):
if stagepos_x[i] > (max(stagepos_x) - (0.3 * xrange)) and turningzone is False:
turncount = turncount + 1
turningzone = True
elif stagepos_x[i] < (max(stagepos_x) - (0.3 * xrange)) and turningzone is True:
turningzone = False
elif stagepos_x[i] < (min(stagepos_x) + (0.3 * xrange)) and turningzone is False:
turncount = turncount + 1
turningzone = True
elif stagepos_x[i] > (min(stagepos_x) + (0.3 * xrange)) and turningzone is True:
turningzone = False
log_print(f'Number of lines: {turncount+1}')
[docs]
class Disk:
def __init__(self, x, y, diameter, frame_index):
self.x = x
self.y = y
self.radius = diameter / 2
self.frame_index = frame_index # New attribute for frame index
def __str__(self):
return f"Disk(frame_index={self.frame_index}, x={self.x}, y={self.y}, radius={self.radius})"
[docs]
def distance(disk1, disk2):
return math.sqrt((disk1.x - disk2.x)**2 + (disk1.y - disk2.y)**2)
[docs]
def are_touching(disk1, disk2):
return distance(disk1, disk2) <= disk1.radius + disk2.radius
[docs]
def find_clusters(disks):
graph = defaultdict(list)
n = len(disks)
for i in range(n):
for j in range(i+1, n):
if are_touching(disks[i], disks[j]):
graph[i].append(j)
graph[j].append(i)
def dfs(node, visited, component):
visited[node] = True
component.append(node)
for neighbor in graph[node]:
if not visited[neighbor]:
dfs(neighbor, visited, component)
visited = [False] * n
clusters = []
for i in range(n):
if not visited[i]:
component = []
dfs(i, visited, component)
clusters.append(component)
return clusters
[docs]
def plot_clusters(input_path, disks, clusters):
fig, ax = plt.subplots()
# Find strongest frame in each cluster
seeds = find_seeds(input_path, clusters)
# Initialize variables to determine the plot limits
min_x, max_x = float('inf'), float('-inf')
min_y, max_y = float('inf'), float('-inf')
# Generate a random color for each cluster
cluster_colors = [np.random.rand(3,) for _ in clusters]
for cluster_id, cluster in enumerate(clusters):
# Calculate the mean position of the cluster to place the cluster index label
cluster_x = [disks[disk_index].x for disk_index in cluster]
cluster_y = [disks[disk_index].y for disk_index in cluster]
mean_x = np.mean(cluster_x)
mean_y = np.mean(cluster_y)
# Update plot limits
min_x, max_x = min(min_x, min(cluster_x)), max(max_x, max(cluster_x))
min_y, max_y = min(min_y, min(cluster_y)), max(max_y, max(cluster_y))
# Plot each disk in the cluster with an individual color
for disk_index in cluster:
disk = disks[disk_index]
color = cluster_colors[cluster_id]
circle = plt.Circle((disk.x, disk.y), disk.radius, fill=True, color=color, alpha=1)
ax.add_artist(circle)
# Highlight the seed frame
#if disk_index == seeds[cluster_id]:
# # Draw a larger circle or change color to highlight the seed frame
# highlight_circle = plt.Circle((disk.x, disk.y), disk.radius * 1.2, fill=False, color='red', linewidth=3)
# ax.add_artist(highlight_circle)
# Plot the cluster index as text near the mean position of the cluster
#ax.text(mean_x, mean_y, f'Cluster {cluster_id}', color='black', fontsize=12, ha='center', va='center')
# Adjust the plot limits with some padding
padding = 0.7
x_range, y_range = max_x - min_x, max_y - min_y
ax.set_xlim(min_x - padding * x_range, max_x + padding * x_range)
ax.set_ylim(min_y - padding * y_range, max_y + padding * y_range)
# Set labels and title
ax.set_xlabel('X Coordinate')
ax.set_ylabel('Y Coordinate')
plt.title('Clusters of Data Points')
plt.axis('equal')
plt.show()
[docs]
def write_cluster_ids(input_path, clusters, disks):
with h5py.File(input_path, 'r') as file:
nPeaks = file['entry/data/nPeaks'][:]
cluster_ids = [int(-1)] * len(nPeaks)
for cluster_id, cluster in enumerate(clusters):
for disk_index in cluster:
frame_index = disks[disk_index].frame_index # Access the frame index from the Disk object
cluster_ids[frame_index] = cluster_id
cluster_ids_np = np.array(cluster_ids)
with h5py.File(input_path, 'a') as file: # Open the file in append mode
# Create a dataset for cluster_ids or overwrite if it already exists
if 'entry/data/cluster_ids' in file:
del file['entry/data/cluster_ids'] # Delete the old dataset
file.create_dataset('entry/data/cluster_ids', data=cluster_ids_np)
[docs]
def find_touching_disks(input_path, threshold, beam_diameter, interlaceoffset = None):
with h5py.File(input_path, 'r') as workingfile:
stagepos_x = workingfile['entry/data/stagepos_x_refined'][:]
stagepos_y = workingfile['entry/data/stagepos_y'][:]
nPeaks = workingfile['entry/data/nPeaks'][:]
streakdirection = workingfile['entry/data/streakdirection'][:]
if interlaceoffset is not None:
stagepos_y[streakdirection == -1] += interlaceoffset
disks = []
for i, (x, y, n) in enumerate(zip(stagepos_x, stagepos_y, nPeaks)):
if n > threshold:
disks.append(Disk(x, y, beam_diameter, i))
clusters = find_clusters(disks)
return disks, clusters
[docs]
def find_seeds(input_path, clusters):
# Find the strongest frame in each cluster (=seed), use it later to grow the cristal starting from the seed
with h5py.File(input_path, 'r') as file:
nPeaks = file['entry/data/nPeaks'][:]
seeds = []
for cluster in clusters:
max_peaks = -1
seed_frame = None
for frame in cluster:
if nPeaks[frame] > max_peaks:
max_peaks = nPeaks[frame]
seed_frame = frame
seeds.append(seed_frame)
return seeds
[docs]
def interlacecorrection(h5file_path, zdim, offset = None, on_pick_event_handler=None):
with h5py.File(h5file_path, 'r') as workingfile:
zpath = 'entry/data/' + str(zdim)
stagepos_x = workingfile['entry/data/stagepos_x_refined'][:]
stagepos_y = workingfile['entry/data/stagepos_y'][:]
zdimension = workingfile[zpath][:]
streakdirection = workingfile['entry/data/streakdirection'][:]
# Apply offset to stagepos_x where streakdirection is -1
stagepos_x[streakdirection == -1] += offset
# Determine the scaling and translation factors for the coordinates
x_min, x_max = stagepos_x.min(), stagepos_x.max()
y_min, y_max = stagepos_y.min(), stagepos_y.max()
x_range = x_max - x_min
#print('x range: ' + str(x_max - x_min))
#print(x_min, x_max)
#print('y range: ' + str(y_max - y_min))
#print(y_min, y_max)
# Normalize and visualize
plt.figure()
plt.scatter(stagepos_x, stagepos_y, c=zdimension, cmap='inferno_r', marker='.', s=100)
#plt.colorbar(label='Number of Peaks')
plt.title('Reconstructed Image')
plt.axis('equal') # Ensure equal scaling of the x and y axes
plt.savefig('intensity_map.png')
#plt.show()
scatter = plt.scatter(stagepos_x, stagepos_y, c=zdimension, cmap='inferno_r', marker='.', s=200, picker=True)
plt.colorbar(scatter, label='Number of Peaks')
if on_pick_event_handler is not None:
# Connect the pick event handler
plt.gcf().canvas.mpl_connect('pick_event', on_pick_event_handler)
plt.show()
[docs]
def get_line_bounds(input_path, frame_index, zdim="nPeaks"):
with h5py.File(input_path, 'r') as workingfile:
streakdirection = workingfile['entry/data/streakdirection'][:]
# Validate frame_index
if frame_index >= len(streakdirection) or frame_index < 0:
raise ValueError("frame_index is out of bounds")
# Find the beginning and end of the line in the specified frame
current_direction = streakdirection[frame_index]
start_index = frame_index
while start_index > 0 and streakdirection[start_index - 1] == current_direction:
start_index -= 1
end_index = frame_index
while end_index < len(streakdirection) - 1 and streakdirection[end_index + 1] == current_direction:
end_index += 1
# Find the beginning and end of the previous line, if exists
prev_line_start, prev_line_end = None, None
if start_index > 0:
prev_line_direction = streakdirection[start_index - 1]
prev_line_end = start_index - 1
prev_line_start = prev_line_end
while prev_line_start > 0 and streakdirection[prev_line_start - 1] == prev_line_direction:
prev_line_start -= 1
# Find the beginning and end of the next line
next_line_start, next_line_end = None, None
if end_index < len(streakdirection) - 1:
next_line_direction = streakdirection[end_index + 1]
next_line_start = end_index + 1
next_line_end = next_line_start
while next_line_end < len(streakdirection) - 1 and streakdirection[next_line_end + 1] == next_line_direction:
next_line_end += 1
return (prev_line_start, prev_line_end), (start_index, end_index), (next_line_start, next_line_end)
# current_line_bounds, next_line_bounds = get_line_bounds(input_path)
# print("Current line bounds:", current_line_bounds)
# print("Next line bounds:", next_line_bounds)
[docs]
def has_common_angle(seed_frame_index, neighbor_frame_index, hdf5_data, angle_tolerance=5):
# Retrieve the top 5 angles for the seed and neighbor frames
seed_angles = hdf5_data['entry/data/index_angles'][seed_frame_index, :5]
neighbor_angles = hdf5_data['entry/data/index_angles'][neighbor_frame_index, :5]
# Check for common angles within the tolerance
for angle in seed_angles:
if angle != -1 and any(abs(angle - other_angle) <= angle_tolerance for other_angle in neighbor_angles if other_angle != -1):
return True
return False
[docs]
def refine_clusters(input_path, clusters, cutoff, angle_tolerance):
with h5py.File(input_path, 'r') as file:
nPeaks = file['entry/data/nPeaks'][:]
index_angles = file['entry/data/index_angles'][:]
refined_clusters = []
for cluster in clusters:
# Find the seed frame (strongest frame) in each cluster
seed_frame = max(cluster, key=lambda frame: nPeaks[frame])
seed_angles = index_angles[seed_frame, :5]
# Filter the cluster
refined_cluster = [seed_frame] # Start the cluster with the seed frame
to_process = set(cluster) - {seed_frame} # Remaining frames to process
while to_process:
current_frame = to_process.pop()
current_angles = index_angles[current_frame, :5]
# Check if the current frame shares an angle with any frame in the refined cluster
if any(angle != -1 and any(abs(angle - other_angle) <= angle_tolerance for other_angle in current_angles if other_angle != -1) for frame in refined_cluster for angle in index_angles[frame, :5]):
refined_cluster.append(current_frame)
else:
to_process -= {current_frame} # Remove frame if it doesn't share an angle
if len(refined_cluster) > 1: # Only add non-trivial clusters
refined_clusters.append(refined_cluster)
return refined_clusters