Optimizing High-Throughput SEM for Large-area Defect Characterization in AM Steel
Contents
%matplotlib ipympl
from ipywidgets import Layout, VBox, HBox, Dropdown, Output
from IPython.display import display
import numpy as np
from cv2 import resize, imread, fitEllipse, ellipse, line
import matplotlib.pyplot as plt
import math
import pyarrow.feather as feather
from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
def get_data_plots(df):
pores_df=df[df["Pore Identity"]=="Pore"]
crack_df=df[df["Pore Identity"]=="Crack"]
#bar chart data
bar_chart_data=[[sum(df["Pore Identity"]=="Crack"),sum(df["Pore Identity"]=="Pore")],[round(df[df["Pore Identity"]=="Crack"]["Pore Area"].mean()),round(df[df["Pore Identity"]=="Pore"]["Pore Area"].mean())]]
# Convert area values to log base 10
log_cracks = np.log10(crack_df["Pore Area"])
log_gas_pores = np.log10(pores_df["Pore Area"])
# Calculate histograms
hist_cracks, bins_cracks = np.histogram(log_cracks, bins='auto')
hist_gas_pores, bins_gas_pores = np.histogram(log_gas_pores, bins='auto')
# Create a list of tuples (histogram, bins, label, color)
hist_data = [
(hist_cracks, bins_cracks, 'Cracks', 'blue'),
(hist_gas_pores, bins_gas_pores, 'Pores', 'green'),
]
#rose data
rose_data=[crack_df["Pore Orientation"],pores_df["Pore Orientation"]]
origin = [1000,1000]
# Calculate distances from the origin
def calculate_distances(coordinates, origin):
return [np.sqrt((x - origin[0]) ** 2 + (y - origin[1]) ** 2) for x, y in coordinates]
distances_cracks = calculate_distances(crack_df["Pore Centroid"], origin)
distances_gas_pores = calculate_distances(pores_df["Pore Centroid"], origin)
#scatter data
scatter=[[distances_cracks,distances_gas_pores],[log_cracks,log_gas_pores]]
#bubble data
df["Category"] = np.where(df["Pore Roundness"] >= 0.3, "Pores", "Cracks")
df["Pore Area"] *= 1
max_area = df["Pore Area"].max()
area=[crack_df["Pore Area"],pores_df["Pore Area"]]
roundness=[crack_df["Pore Roundness"],pores_df["Pore Roundness"]]
bubble=[area,roundness,max_area]
return bar_chart_data,hist_data,rose_data,scatter,bubble
def bar_chart(ax,data):
counts,averages=data
colors = ['blue', 'green']
labels = ['Cracks', 'Pores']
bars = ax.bar(labels, counts, color=colors)
# Add the average size labels on the bars
for bar, avg in zip(bars, averages):
yval = bar.get_height()
ax.text(bar.get_x() + bar.get_width()/2, yval+0.05*max(counts), f'Average Size: {avg:.0f}', ha='center', va='top', color='black', fontsize=7)
# Add labels and title
ax.set_ylabel('Number of Defects')
ax.set_title('Defects Amount and Their Average Size', fontsize=9)
def histogram(ax,hist_data):
# Sort the histograms by the maximum frequency
hist_data.sort(key=lambda x: max(x[0]), reverse=True)
# Plot histograms
for hist, bins, label, color in hist_data:
ax.hist(bins[:-1], bins, weights=hist, alpha=0.7, label=label, color=color)
# Add labels and title
ax.set_xlabel('Log base 10 of Area')
ax.set_ylabel('Frequency')
ax.set_title('Area for Different Defect Types', fontsize=10)
ax.legend()
def rose_diagram(axs,data,cracks):
def helper(ax, data, title, color, bin_edges, agreements=None, annotate=False):
angles_radians = np.deg2rad(data)
hist, _ = np.histogram(angles_radians, bins=bin_edges)
bin_widths = np.diff(bin_edges)
bars = ax.bar(bin_edges[:-1], hist, width=bin_widths, edgecolor='black', color=color, alpha=0.6)
ax.set_theta_zero_location('E')
ax.set_theta_direction(1)
ax.set_thetalim(0, np.pi)
ax.set_xticks(np.linspace(0, np.pi, len(bin_edges) // 2 + 1))
ax.set_xticklabels([f'{int(np.rad2deg(angle))}°' for angle in np.linspace(0, np.pi, len(bin_edges) // 2 + 1)])
ax.set_title(title)
# Reduce number of radial ticks (e.g., only show 0, 5, and 10)
ax.set_yticks([0, 50, 100, 140])
# Optionally, reduce the font size of radial labels
ax.tick_params(axis='y', labelsize=8)
if annotate and agreements is not None:
overall_agreement_5to10, overall_agreement_1to10, bin_agreement_5to10, bin_agreement_1to10 = agreements
ax.text(0.5, 1.1, f'Agreement 5ys-10ys: {overall_agreement_5to10:.2f}\nAgreement 1ys-10ys: {overall_agreement_1to10:.2f}', transform=ax.transAxes, ha='center', fontsize=10)
for i, (agreement_5to10, agreement_1to10) in enumerate(zip(bin_agreement_5to10, bin_agreement_1to10)):
angle = bin_edges[i] + (bin_edges[i+1] - bin_edges[i]) / 2
ax.text(angle, max(hist) * 0.7, f'{agreement_5to10:.2f}\n{agreement_1to10:.2f}', color='blue', ha='center', fontsize=11)
num_bins = 18
bin_edges = np.linspace(0, np.pi, num_bins + 1)
if cracks:
helper(axs,data[0],"Cracks","blue",bin_edges)
else:
helper(axs,data[1],"Pores","red",bin_edges)
def scatter_plot(ax,data):
distances,log=data
ax.scatter(distances[0], log[0], color='blue', alpha=0.6, edgecolors="w", linewidth=0.5, label='Cracks')
ax.scatter(distances[1], log[1], color='green', alpha=0.6, edgecolors="w", linewidth=0.5, label='Pores')
ax.set_title('Area vs Distance from Origin')
ax.set_xlabel('Distance from Origin')
ax.set_ylabel('Log Area')
ax.legend()
ax.grid(True)
def bubble_plot(axs,data):
area,roundness,max=data
bubble_size_scale = 700
colors = {"Cracks": "blue", "Pores": "green"}
i=0
for category, color in colors.items():
axs.scatter(
area[i],
roundness[i],
s=(np.sqrt(area[i]) / np.sqrt(max)) * bubble_size_scale,
c=color,
alpha=0.7,
edgecolors="w",
linewidth=0.5,
label=category
)
i+=1
axs.set_xlabel("Pore Area (scaled μm²)")
axs.set_ylabel("Roundness")
axs.set_title("Bubble Plot: Area vs Roundness", fontsize=9)
axs.grid(True)
axs.legend(title="Categories")
def on_click(event,axs,filled_detection,mini_axs):
message=(
"not a pore "
)
rows,cols=filled_detection.shape
if event.inaxes == axs:
j = round(event.xdata)
i = round(event.ydata)
if filled_detection[i, j] <=0:
filled_detection[i, j] = 10
Marked_List = []
Unmarked_List=[[i,j]]
List_Boundary=[]
P = 0
while len(Unmarked_List)>0:
k1,k2 = Unmarked_List[-1]
Marked_List.append(Unmarked_List[-1])
Unmarked_List.pop(-1)
for dx, dy in [(-1, 0), (1, 0), (0, -1), (0, 1)]:
nx, ny = k1 + dx, k2 + dy
if 0 <= nx < rows and 0 <= ny < cols and filled_detection[nx, ny] == 0:
filled_detection[nx, ny] = 10
Unmarked_List.append([nx, ny])
elif (0 <= nx < rows and 0 <= ny < cols)==False:
P = P + 1
List_Boundary.append([k1,k2])
elif filled_detection[nx, ny] == 200:
P = P + 1
List_Boundary.append([k1,k2])
A = len(Marked_List)
for x,y in Marked_List:
filled_detection[x,y]=0
try:
List_Boundary=np.array(List_Boundary)
roundness = 4 * math.pi * A / (P**2)
if len(List_Boundary)<5:
pass
elif roundness>0.3 :#ellipse fiiting 0.3
fitted_ellipse = fitEllipse(List_Boundary)
center, axes, angle = fitted_ellipse
fitted_ellipse=((center[1],center[0]),axes,angle)
angle=(90-angle)
if angle<0:
angle+=180
message = (
f"📌 **Defect Details** ⚪ Roundness: {roundness:.3f}\n"
f"📍 Centroid: {int(center[0]),int(center[1])} 📐 Orientation: {angle}°\n"
f"📏 Size: {A:.2f} μm² 🔹 Classification: Ellipse"
)
separation=round((A/3)**0.5)+10
mini_axs.clear()
filled_detection3=np.copy(filled_detection)
ellipse(filled_detection3, fitted_ellipse, (255, 255, 255), 10)
for x,y in Marked_List:
filled_detection3[x,y]=10
mini_axs.imshow(filled_detection3[i-separation:i+separation,j-separation:j+separation])
plt.draw()
else:
#skeleton
# Convert coordinates to numpy array
coordinates_array = np.array(List_Boundary)
# Calculate pairwise distances between coordinates
distances = np.linalg.norm(coordinates_array[:, None] - coordinates_array, axis=-1)
# Find the indices of the two points with the maximum distance
idx_max = np.unravel_index(np.argmax(distances), distances.shape)
point1, point2 = List_Boundary[idx_max[0]], List_Boundary[idx_max[1]]
angle=round(math.atan(float(-(point1[0]-point2[0])/(point1[1]-point2[1])))*180/math.pi)
if angle<0:
angle+=180
angle=angle%180
message = (
f"📌 **Defect Details** ⚪ Roundness: {roundness:.3f}\n"
f"📍 Centroid: {(int(point1[0]+point2[0])/2),(int(point1[1]+point2[1])/2)} 📐 Orientation: {int(angle)}°\n"
f"📏 Size: {A:.2f} μm² 🔹 Classification: Crack"
)
height=round(abs(point1[0]-point2[0])/2)+10
width=round(abs(point1[1]-point2[1])/2)+10
mini_axs.clear()
filled_detection3=np.copy(filled_detection)
line(filled_detection3, (point2[1],point2[0]), (point1[1],point1[0]),(255, 255, 255),1)
mini_axs.imshow(filled_detection3[i-height:i+height,j-width:j+width])
plt.draw()
except ZeroDivisionError:
pass
with output_text:
output_text.clear_output(wait=True)
print(message)
class DataAnalysisWidget:
def __init__(self, image):
self.image = image # Imagen cargada
self.out = Output(layout=Layout()) # Salida para los gráficos display="flex", justify_content='flex-end', margin='0 0 0 80px'
self.mini_out=Output(layout=Layout())
with self.out:
self.fig, self.axs = plt.subplots(figsize=(3, 3))
self.polar_axs=None
self.axs.imshow(image)
self.handler=self.fig.canvas.mpl_connect(
"button_press_event", lambda event: on_click(event, self.axs, image,self.mini_axs)
)
scalebar = AnchoredSizeBar(self.axs.transData,
250, '50 μm', 'lower right',
pad=0.1,
color='black',
frameon=False,
size_vertical=1)
self.axs.add_artist(scalebar)
self.fig.set_label(' ')
self.axs.set_title("Analysed Image")
self.axs.set_aspect("auto") # Adjust automatically
plt.show()
with self.mini_out:
self.mini_fig, self.mini_axs = plt.subplots(figsize=(3, 3))
self.mini_axs.imshow(np.zeros((500,500)))
self.mini_fig.canvas.toolbar_visible=False
self.mini_fig.canvas.callbacks.callbacks = {} # Remove all event handlers
self.mini_fig.tight_layout()
self.mini_fig.set_label(' ')
self.mini_axs.set_title("Orientation Plot")
self.mini_fig.subplots_adjust(left=0.2, right=0.95, bottom=0.2, top=0.9)
plt.show()
def update_plot(self,event):
value = event["new"]
self.fig.canvas.mpl_disconnect(self.handler)
self.mini_axs.clear()
self.mini_axs.set_aspect("auto")
message=""
if self.polar_axs is not None:
self.mini_fig.delaxes(self.polar_axs)
self.polar_axs=None
if value == "Scatter":
# Generar un gráfico de dispersión
scatter_plot( self.mini_axs,scatter)
message= """
Displays defects where the
x-axis shows distance from
a reference point and the
y-axis represents the logarithm
of defect area. Blue points
indicate cracks and green points
indicate pores.
"""
elif value == "Bar Chart":
bar_chart( self.mini_axs,bar_chart_data)
message= """
Shows the total number of defects
along with their average sizes,
with blue bars for cracks and
green bars for pores.
"""
elif value=="Pores Rose":
self.mini_axs.set_axis_off()
self.polar_axs=self.mini_fig.add_subplot(111, projection='polar')
rose_diagram(self.polar_axs,rose_data,cracks=False)
message= """
A polar plot that maps the
orientation distribution of
pores using angular bins (angle
measured from the horizontal).
"""
elif value == "Bubble":
# Gráfico de burbujas (intensidad como tamaño)
bubble_plot( self.mini_axs,bubble)
message= """
Correlates defect area with
roundness; bubble sizes are scaled
by area, with blue for cracks and
green for pores
"""
elif value == "Histogram":
# Histograma de intensidades
histogram( self.mini_axs,hist_data)
message= """
Displays the frequency
distribution of defect areas on
a logarithmic scale, differentiating
between cracks (blue) and
pores (green).
"""
elif value == "General":
# Mostrar la imagen original
self.axs.imshow(self.image)
self.mini_axs.imshow(np.zeros((500,500)))
self.axs.set_title("General Plot")
self.handler=self.fig.canvas.mpl_connect(
"button_press_event", lambda event: on_click(event, self.axs, self.image,self.mini_axs)
)
with output_text:
output_text.clear_output(wait=True)
print(message)
plt.draw()
output_text=Output(layout=Layout(height="100px"))
path = "../Data/10ys Defects.npz"
img = np.load(path)["arr_0"][:,:]
path_data="../Data/DataforWidget2.feather"
df = feather.read_feather(path_data)
bar_chart_data,hist_data,rose_data,scatter,bubble=get_data_plots(df)
# Crear lista de opciones para el dropdown
plot_list = ["General", "Histogram", "Bubble","Pores Rose", "Bar Chart", "Scatter"]
dropdown_widget = Dropdown(
options=plot_list, value="General", description="Plots", layout=Layout(width="200px")
)
d=DataAnalysisWidget(img)
dropdown_widget.observe(lambda event: d.update_plot(event), names="value")
v=VBox([output_text,dropdown_widget],layout=Layout(align_items='center'))
h=HBox([d.out,d.mini_out])
V=VBox([h,v], layout=Layout(align_items='center'))
display(V)
VBox(children=(HBox(children=(Output(), Output())), VBox(children=(Output(layout=Layout(height='100px')), Drop…
