Redesigned a couple function to better use the functionality of pytho…

CD_SEM_Calc.py

@@ -2,6 +2,8 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.fftpack import fftshift, fft2  # , ifft2
+from scipy.ndimage import median_filter
+from skimage.draw import line
 
 # from scipy.ndimage import gaussian_filter, morphology
 from scipy.stats import scoreatpercentile
@@ -29,7 +31,7 @@ def image_size(
     imax = int(2 ** np.floor(np.log2(height)))
     lmax = (imax + 40) if height >= (imax + 40) else imax
     kscale = 2 * np.pi / (lmax * pixel_size * pixel_scale * 10**3)
-    return imax, lmax, kscale
+    return imax, (lmax-1)|1, kscale # Guarantees lmax is ODD insuring a pixel at the very center of the image
 
 
 def extract_center_part(img: np.ndarray, size: int) -> np.ndarray:
@@ -44,13 +46,14 @@ def extract_center_part(img: np.ndarray, size: int) -> np.ndarray:
         np.arry: subarray defined by the indices
     """
     height, width = img.shape
-    roi_h = int(np.maximum(np.floor((height - size) / 2), 0))
-    roi_w = int(np.maximum(np.floor((width - size) / 2), 0))
+    dimension = min(height, width)
+    roi_h = int(max(0,np.floor((dimension - size) / 2)))
+    roi_w = int(max(0,np.floor((dimension - size) / 2)))
     # roi = 0  when there is a data zone below zero
     return img[roi_h : int(roi_h + size), roi_w : int(roi_w + size)]
 
 
-def fourier_img(img: np.ndarray) -> tuple[np.ndarray, float]:
+def fourier_img(img: np.ndarray) -> tuple[np.ndarray, tuple]:
     """Calculates the magnitude squared of the Fourier Transform of a square image "img" of size "lmax".
     It recenters it so that the zero frequency is at {lmax/2+1, lmax/2+1}. It saves the magnitude square
     of the zero frequency in a variable called "ctrval". The zero frequncy in the image is replaced by "1" to help
@@ -60,17 +63,20 @@ def fourier_img(img: np.ndarray) -> tuple[np.ndarray, float]:
         img (np.ndarray): Clipped and Rescaled image that need processed
 
     Returns:
-        tuple[np.ndarray, float]: FFT image, Magnitude square of the zero frequency
+        tuple[np.ndarray, tuple]: FFT image, Magnitude square of the zero frequency
     """
-    fimg = np.abs(fftshift(fft2(img))) ** 2
+    fimg = np.abs(fftshift(fft2(np.copy(img)))) ** 2
     center = np.array(fimg.shape) // 2
-    ctrval = fimg[center, center]
-    fimg[center, center] = 1
-    return tools.rescale_array(np.log(fimg), 0, 1), ctrval
+    rescale_values = (fimg[center, center], np.min(fimg), np.max(fimg))
+    fimg = tools.rescale_array(np.log(fimg), 0, 1)
+    return fimg, rescale_values
 
 
-def rotated_angle(probe: int, img: np.ndarray, lmax: int) -> float:
-    """_summary_
+def rotated_angle(angle: float, img: np.ndarray, lmax: int) -> float:
+    """Integrates the I(qx, qy = 0) line (the horizontal center line).
+    It also probes different rotation values by rotating the image one pixel at a time up to "probe" pixels.
+    It then plots the integral values vs rotated amount. It picks the maximum value to be the one that indicates the "true" horizontal position.
+    It returns the rotated angle in degrees. "probe needs to be an ODD number"
 
     Args:
         probe (int): Maximum number of pixels to rotate by
@@ -80,38 +86,21 @@ def rotated_angle(probe: int, img: np.ndarray, lmax: int) -> float:
     Returns:
         float: angle the image needs rotated
     """
+    def line_sum(a, img):
+        dy = int(np.ceil((lmax/2)*np.arctan(np.radians(a))))
+        rr, cc = line(0, int(lmax//2 - dy), lmax-1, int(lmax//2 + dy))
+        return np.sum(img[rr, cc])
 
-    def calculatetotals(probe: int, img: np.ndarray, lmax: int) -> np.ndarray:
-        totals = []
-        for l in range(-probe, probe + 1):
-            indices = np.array(
-                [
-                    (round((j - lmax / 2 - 1) * (l / (lmax / 2 - 1)) + lmax / 2 + 1), j)
-                    for j in range(lmax)
-                ]
-            )
-            values = img[indices[:, 0], indices[:, 1]]
-            total = np.sum(values)
-            totals.append(total)
-        totals = np.array(totals)
-        return totals
-
-    def movingmedian(data: np.ndarray, window_size: int) -> float:
-        padded_data = np.pad(data, (window_size - 1) // 2, mode="edge")
-        windowed_data = np.lib.stride_tricks.sliding_window_view(
-            padded_data, window_size
-        )
-        medians = np.apply_along_axis(lambda x: np.median(x), 1, windowed_data)
-        return medians
-
-    totals = calculatetotals(probe, img, lmax)
-    totals2 = movingmedian(totals, 7)
+    probe = int(np.ceil((lmax/2)*np.arctan(np.radians(angle))))
+    totals = [line_sum(n, img) for n in range(-angle, angle + 1)]
+    totals2 = median_filter(totals, size=3, mode='mirror')
 
     max_index = np.argmax(totals2)
-    if totals2[max_index] / totals[probe] > 1.05:
-        ra = np.arctan((max_index - (probe - 3 + 1)) / (lmax / 2 - 1)) * 180 / np.pi
-    else:
-        ra = 0
+    ra =0
+    # if totals2[max_index] / totals[probe] > 1.05:
+    #     ra = np.arctan((max_index - (probe - 3 + 1)) / (lmax / 2 - 1)) * 180 / np.pi
+    # else:
+    #     ra = 0
     print("Angle of rotation:", ra)
 
     angle_range = np.linspace(-probe, probe, len(totals))

CD_SEM_tools.py

@@ -29,3 +29,19 @@ def rescale_array(arr: np.ndarray, new_min: float, new_max: float) -> np.ndarray
         scaled_arr * (new_max - new_min) + new_min
     )  # Rescale to the new range
     return rescaled_arr
+
+
+def threashold_mask(img: np.ndarray, threshold: float, new_value: float = 0) -> np.ndarray:
+    """Sets all values in an array below the threashold to the new_value
+
+    Args:
+        img (np.ndarray): imgage array to have the mask applied to
+        threshold (float): The lowest value being kept
+        new_value (float): The value assigned to all pixels below threshold: Default zero
+
+    Returns:
+        np.ndarray: Output image array
+    """
+    new_img = np.copy(img)
+    new_img[img < threshold] = new_value
+    return new_img

SEM_Image.py

@@ -11,6 +11,7 @@
 TAG_NAMES: Final[list[str]] = ["ImageWidth", "ImageLength", "CZ_SEM"]
 PIX_SIZE: Final[str] = "ap_image_pixel_size"
 IMAGE_SCALE_UM: Final[float] = 0.2  # Image scale bar length in micrometers
+MASK_THRESHOLD: Final[float] = 0.9  # Threshold for binary mask
 
 ######### This object holds 54 properties from CD-SEM analyis. It only auto initializes values pulled directly from the header file of the '.fit' SEM image.
 ######### The __call__ function will run all the necassary calculations to assign values to all object properties
@@ -27,12 +28,8 @@ def __init__(self):
         self.imax: int | None = None  # See ImgFlat.lmax for details on variables
         self.lmax: int | None = None
         self.kscale: float | None = None
-        self.image_FFT_center: None | float = (
-            None  # Magnitude square of the zero frequency of the FFT image
-        )
-        self.image_rotate: None | float = (
-            None  # The angle the image needs rotated to make sure the FFT is horizontal
-        )
+        self.rescale_values: None | list = [None] * 3  # Magnitude square of the zero frequency, minimum and maximum pixel values from the unnormalized FFT image
+        self.image_rotate: None | float = None  # The angle the image needs rotated to make sure the FFT is horizontal
 
         # Images
         self.image = tifffile.imread(self.path)  # Original
@@ -91,14 +88,13 @@ def __init__(self):
         self.BLPA_inline: None | float = None  # nm # In Line
 
     def __call__(self):
+
         self.imax, self.lmax, self.kscale = calc.image_size(
             self.height, self.pix_scale, self.pix_size
         )
-        self.image = tools.rescale_array(
-            calc.extract_center_part(self.image, self.lmax), 0, 1
-        )
-        self.image_FFT, self.image_FFT_center = calc.fourier_img(self.image)
-        self.image_rotate = calc.rotated_angle(25, self.image_FFT, self.lmax)
+        self.image = calc.clip_image(calc.extract_center_part(self.image, self.lmax))
+        self.image_FFT, self.rescale_values = calc.fourier_img(np.copy(self.image))
+        self.image_rotate = calc.rotated_angle(25, tools.threashold_mask(self.image_FFT, MASK_THRESHOLD), self.lmax)
 
     def _sem_image_selector(self) -> str:
         """Lets you select the image file for the object
@@ -172,12 +168,11 @@ def display_SEM_image(self: object, image: tifffile, bar=False) -> None:
 
         if bar:
             # Calculate the dimensions of the scale bar
-            image_height = image.shape[0]
-            scale_bar_length_pixels = (IMAGE_SCALE_UM * self.pix_scale) / self.pix_size
+            scale_bar_length_pixels = IMAGE_SCALE_UM / (self.pix_scale * self.pix_size)
 
             # Calculate the position of the scale bar
             scale_bar_x = image.shape[1] - scale_bar_length_pixels - 100
-            scale_bar_y = image_height - 50
+            scale_bar_y = image.shape[0] - 50
 
             # Add the scale bar to the plot
             scale_bar = patches.Rectangle(
@@ -204,19 +199,25 @@ def display_SEM_image(self: object, image: tifffile, bar=False) -> None:
         # Show the plot
         plt.show()
 
-    def display_fft_image(self: object, fimg: np.ndarray) -> None:
+    def display_fft_image(self: object, fimg: np.ndarray, mask=False) -> None:
         """Displays the scaled FFT image on a colorblind friendly colorbar
 
         Args:
             fimg (np.ndarray): FFT image getting displayed
         """
         # Define a colorblind-friendly colormap
+        global MASK_THRESHOLD
+
         cmap = LinearSegmentedColormap.from_list(
             "colorblind_cmap", ["#000000", "#377eb8", "#ff7f00", "#4daf4a"], N=256
         )
 
+        if mask:
+            plt.imshow(tools.threashold_mask(fimg, MASK_THRESHOLD), cmap=cmap)
+        else:
+            plt.imshow(fimg, cmap=cmap)
+
         # Plot the FFT image
-        plt.imshow(fimg, cmap=cmap)
         plt.colorbar(label="Intensity")
         plt.title("FFT Image")
         plt.axis("off")