Optimizing High-Throughput SEM for Large-area Defect Characterization in AM Steel

%matplotlib ipympl

from ipywidgets import Layout, VBox, Output, Button, FloatSlider
from IPython.display import display
import numpy as np
from cv2 import adaptiveThreshold, ADAPTIVE_THRESH_GAUSSIAN_C,THRESH_BINARY, Canny, dilate, getStructuringElement, MORPH_ELLIPSE
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import math
from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
from scipy import ndimage as ndi
from scipy.ndimage import gaussian_filter


image_quality_dpi=300

# Widget 1 - takes the Raw Image and Displays the Gaussian Detection in Binary #

def Widget_1_Function(Raw_Image,fidelity=40,max_range=151):
    Gaussian_Detection = np.full_like(Raw_Image, 200)
    gaussian_ranges = [31, max_range]

    for gaussian_range in gaussian_ranges:
        new_raw = Raw_Image.copy()
        gaussian_filter = adaptiveThreshold(Raw_Image, 200, ADAPTIVE_THRESH_GAUSSIAN_C, THRESH_BINARY,
                                               gaussian_range, fidelity)
        Gaussian_Detection[gaussian_filter == 0] = 0
        new_raw[gaussian_filter == 0] = 255
        del new_raw

    return(Gaussian_Detection)

# Widget 2 - takes the Raw Image and Displays the Canny Detection #

def Widget_2_Function(Raw_Image,low_t=11,high_t=70,sigma=3):
    """Recalculate edges whenever any slider is moved."""

    # Apply Gaussian filter with current sigma
    blurred = gaussian_filter(Raw_Image, sigma=sigma)
    blurred_uint8 = np.uint8(blurred)

    # Perform Canny
    Canny_Detection = Canny(blurred_uint8, low_t, high_t)
    Canny_Detection[Canny_Detection == 0] = 200
    Canny_Detection[Canny_Detection == 255] = 0
    return(Canny_Detection)

# Widget 3 - takes the Combined Detection (Gaussian + Canny), opens it and then fills #

def Widget_3_Function(Canny_Detection,Gaussian_Detection):
    combined_detection=np.zeros(np.shape(Canny_Detection))
    combined_detection[combined_detection==0]=200
    combined_detection[Canny_Detection == 0] = 0
    combined_detection[Gaussian_Detection == 0] = 0
    inverted_array = np.where(combined_detection == 200, 0, np.where(combined_detection == 0, 200, combined_detection))
    filled_detection = ndi.binary_fill_holes(inverted_array).astype(int)
    filled_detection = np.where(filled_detection == 200, 0,
                                np.where(filled_detection == 0, 200, filled_detection)).astype(np.uint8)
    filled_detection = dilate(filled_detection, getStructuringElement(MORPH_ELLIPSE, (5, 5)), iterations=1)
    filled_detection[Gaussian_Detection == 0] = 0
    filled_detection[filled_detection == 1] = 0

    inverted_array = np.where(filled_detection == 200, 0, np.where(filled_detection == 0, 200, filled_detection))
    filled_detection = ndi.binary_fill_holes(inverted_array).astype(int)
    filled_detection = np.where(filled_detection == 200, 0,
                                np.where(filled_detection == 0, 200, filled_detection)).astype(np.uint8)
    filled_detection[filled_detection == 1] = 0
    Filled_Detection=filled_detection
    return(Filled_Detection)

# Widget 4 - Filters the Filled Detection and Overlaps Results onto the Raw Image #

def Widget_4_Function(Filled_Detection,Raw_Image):
    Contrast_Image=Raw_Image
    Contrast_Image[Contrast_Image>170]=np.mean(Raw_Image)

    pore_cut_off=20

    filled_detection=Filled_Detection
    raw_image=Raw_Image

    final_pores = np.full_like(raw_image, 200)
    rows, cols = raw_image.shape

    raw_image = original_image.copy()

    rows, cols = raw_image.shape

    for i in range(rows):
        for j in range(cols):
            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)

                X_vals = [element[0] for element in Marked_List]

                Y_vals = [element[1] for element in Marked_List]

                try:
                    roundness = 4 * math.pi * A / (P ** 2)
                    pore_value = roundness * A
                    centroid_point = [sum(X_vals) // A, sum(Y_vals) // A]

                    # Calculate distances from boundary points to centroid
                    distances = [np.sqrt((x - centroid_point[0]) ** 2 + (y - centroid_point[1]) ** 2)
                                 for x, y in List_Boundary]

                    # Count extrema in distances
                    extrema_count = 0
                    for k in range(1, len(distances) - 1):
                        if (distances[k] > distances[k - 1] and distances[k] > distances[k + 1]) or \
                                (distances[k] < distances[k - 1] and distances[k] < distances[k + 1]):
                            extrema_count += 1

                    # Calculate the spiky metric
                    spiky_metric = extrema_count / len(distances) if len(distances) > 0 else 0

                    if pore_value >= pore_cut_off and spiky_metric < 0.3 and roundness > 0.05:
                        for x, y in Marked_List:
                            final_pores[x, y] = 0


                except ZeroDivisionError:
                    continue

    Contrast_Image[final_pores==0]=200
    return(Contrast_Image)



path = "../Data/2025 Label for 12-2023 Data - 10ys.npz"
original_image = np.load(path)["arr_0"][:,:,0].astype('uint8')
original_copy=original_image.copy()
Gaussian_Detection=Widget_1_Function(original_copy)
Canny_Detection=Widget_2_Function(original_copy)
Filled_Detection=Widget_3_Function(Gaussian_Detection,Canny_Detection)
Contrast_Image=Widget_4_Function(Filled_Detection,original_copy)

images = [original_image,Gaussian_Detection, Canny_Detection, Filled_Detection, Contrast_Image]
titles=["Raw image","Gaussian Detection","Canny Detection","Filled Detection","Contrast Image"]

# Create the figure
fig = plt.figure()
gs = gridspec.GridSpec(3,4) 

ax1 = fig.add_subplot(gs[1, 0])  # top-left
ax2 = fig.add_subplot(gs[0:2, 1])  # top-center
ax3 = fig.add_subplot(gs[1:3, 1])  # top-right
ax4 = fig.add_subplot(gs[1, 2])  # center (below top-center)
ax5 = fig.add_subplot(gs[1, 3])  # right of top-right


ax=[ax1,ax2,ax3,ax4,ax5]
fig.set_label(" ")
fig.set_figheight(5)
fig.set_figwidth(7)
original_display = ax[0].imshow(images[0], cmap='gray', interpolation='nearest')
gaussian_display = ax[1].imshow(images[1], cmap='gray', interpolation='nearest')
canny_display = ax[2].imshow(images[2], cmap='gray', interpolation='nearest')
filled_display = ax[3].imshow(images[3], cmap='gray', interpolation='nearest')
contrast_display = ax[4].imshow(images[4], cmap='gray', interpolation='nearest')

scalebar = AnchoredSizeBar(ax[0].transData,
                           250, '50 µm', 'lower right', 
                           pad=0.1,
                           color='black',
                           frameon=False,
                           size_vertical=1)

ax[0].add_artist(scalebar)

for axes,title in zip(ax,["Original Image", "Gaussian Detection", "Canny Detection", "Filled Detection", "Contrast Image"]):
    axes.set_title(title, fontsize=10)
#ax.set_xlabel("Length (pixels)")
#ax.set_ylabel("Height (pixels)")
plt.tight_layout()



# Create the sliders
slider1 = FloatSlider(value=25, min=10, max=50, step=5, description='Gaussian Fidelity:', continuous_update=False,style={'description_width': '150px'},layout=Layout(width='400px'))
slider2 = FloatSlider(value=1, min=1, max=100, step=1, description='Canny Low:', continuous_update=False,style={'description_width': '150px'},layout=Layout(width='400px'))
slider3 = FloatSlider(value=1, min=1, max=250, step=1, description='Canny High:', continuous_update=False,style={'description_width': '150px'},layout=Layout(width='400px'))
slider4 = FloatSlider(value=25, min=25, max=175, step=2, description='Gaussian Range:', continuous_update=False,style={'description_width': '150px'},layout=Layout(width='400px'))
slider5 = FloatSlider(value=1, min=1, max=5, step=0.1, description='Canny Sigma:', continuous_update=False,style={'description_width': '150px'},layout=Layout(width='400px'))



def adf(b):
    global images
    original_copy=original_image.copy()
    Gaussian_Detection=Widget_1_Function(original_copy,int(slider1.value),int(round(slider4.value)))
    Canny_Detection=Widget_2_Function(original_copy,int(slider2.value),int(slider3.value),int(slider5.value))
    Filled_Detection=Widget_3_Function(Gaussian_Detection,Canny_Detection)
    Contrast_Image=Widget_4_Function(Filled_Detection,original_copy)

    images = [original_image,Gaussian_Detection, Canny_Detection, Filled_Detection, Contrast_Image]

    original_display.set_data(original_image)
    gaussian_display.set_data(Gaussian_Detection)
    canny_display.set_data(Canny_Detection)
    filled_display.set_data(Filled_Detection)
    contrast_display.set_data(Contrast_Image)


# Create the button
button = Button(description="Update Image")
button.on_click(adf)

# Run the animation
output = Output(layout=Layout(align_items='center'))
right_panel = VBox([slider1, slider4, slider2, slider3, slider5, button], layout=Layout(align_items='center', justify_content='center', background_color="white"))

# Create the main horizontal layout (image on the left, sliders + button on the right)
ui = VBox([VBox([output, right_panel], layout=Layout(align_items='center'))],layout=Layout(align_items="flex-start"))


with output:    
    plt.show()


display(ui)