added evaulation for 2d fft

cal.py

@@ -1,12 +1,12 @@
 import numpy as np
-from scipy.stats import multivariate_normal
 import matplotlib.pyplot as plt
-import scipy
 import scipy.fftpack as sfft
+import matplotlib.patches as patches
+import matplotlib
 
 
 def deg_2_rad(winkel):
-    return winkel/180.0 * np.pi
+    return winkel / 180.0 * np.pi
 
 
 # all units in angstrom
@@ -28,7 +28,6 @@ def mono_2_rutile(c_m, a_m):
 
 
 class Lattice:
-
     def _get_rutile(self, X, Y):
         self.atom_x_rut = X * base_c_r + np.mod(Y, 4) * 0.5 * base_c_r
         self.atom_y_rut = Y * 0.5 * base_a_r
@@ -41,8 +40,7 @@ class Lattice:
 
         print("A_r: ", offset_a_r, "C_r: ", offset_c_r)
 
-        self.atom_x_mono = offset_a_r + X * \
+        self.atom_x_mono = offset_a_r + X * base_c_r + np.mod(Y, 4) * 0.5 * base_c_r
-            base_c_r + np.mod(Y, 4) * 0.5 * base_c_r
         self.atom_x_mono[np.mod(X, 2) == 0] -= 2 * offset_a_r
 
         self.atom_y_mono = offset_c_r + 0.5 * Y * base_a_r
@@ -52,15 +50,15 @@ class Lattice:
         x = np.arange(x_len)
         y = np.arange(y_len)
         X, Y = np.meshgrid(x, y)
-        X[np.mod(Y, 4) == 3] = X[np.mod(Y, 4) == 3]-1
+        X[np.mod(Y, 4) == 3] = X[np.mod(Y, 4) == 3] - 1
-        X[np.mod(Y, 4) == 2] = X[np.mod(Y, 4) == 2]-1
+        X[np.mod(Y, 4) == 2] = X[np.mod(Y, 4) == 2] - 1
         assert np.mod(x.size, 2) == 0
         assert np.mod(y.size, 2) == 0
 
         return X, Y
 
     def __init__(self, x_len: int, y_len: int):
-        X, Y = self._generate_vec(x_len*2, y_len*2)
+        X, Y = self._generate_vec(x_len * 2, y_len * 2)
         self._get_mono(X, Y)
         self._get_rutile(X, Y)
 
@@ -102,6 +100,7 @@ def test_lattice():
 
 
 RESOLUTION = 0.1
+CMAP = "Greys"
 
 
 def image_from_pos(pos_x, pos_y):
@@ -141,40 +140,175 @@ def gaussian(img):
     x = np.linspace(-ratio, ratio, img.shape[0])
     y = np.linspace(-1, 1, img.shape[1])
     X, Y = np.meshgrid(x, y)
-    sigma = 0.3
+    sigma = 0.5
-    z = (1/(2*np.pi*sigma*sigma) * np.exp(-(X**2/(2*sigma**2)
+    z = (
-                                            + Y**2/(2*sigma**2))))
+        1
+        / (2 * np.pi * sigma * sigma)
+        * np.exp(-(X**2 / (2 * sigma**2) + Y**2 / (2 * sigma**2)))
+    )
     return np.multiply(img, z.T)
 
 
+def padding(array, xx, yy):
+    """
+    :param array: numpy array
+    :param xx: desired height
+    :param yy: desirex width
+    :return: padded array
+    """
+
+    h = array.shape[0]
+    w = array.shape[1]
+
+    a = (xx - h) // 2
+    aa = xx - a - h
+
+    b = (yy - w) // 2
+    bb = yy - b - w
+
+    return np.pad(array, pad_width=((a, aa), (b, bb)), mode="constant")
+
+
+def rect_at_point(x, y, color):
+    length_2 = 0.08
+    rect = patches.Rectangle(
+        (x - length_2, y - length_2),
+        2 * length_2,
+        2 * length_2,
+        linewidth=1,
+        edgecolor=color,
+        facecolor="none",
+    )
+    return rect
+
+
+def reci_rutile():
+    x = np.arange(-2, 3)
+    y = np.arange(-2, 3)
+    X, Y = np.meshgrid(x, y)
+    return (X * 0.22 + Y * 0.44).flatten(), (X * 0.349).flatten()
+
+
+def reci_mono():
+    x, y = reci_rutile()
+    return x + 0.1083, y + 0.1719
+
+
+def draw_big_val_rect(img, x, y, x_index, y_index):
+    length_2 = 0.08
+    pos_x_lower = x - length_2
+    pos_x_upper = x + length_2
+
+    pos_y_lower = y - length_2
+    pos_y_upper = y + length_2
+    x_lower = np.searchsorted(x_index, pos_x_lower)
+    x_upper = np.searchsorted(x_index, pos_x_upper)
+    y_lower = np.searchsorted(y_index, pos_y_lower)
+    y_upper = np.searchsorted(y_index, pos_y_upper)
+
+    img[y_lower:y_upper, x_lower:x_upper] = 1e4
+    return img
+
+
+def extract_rect(img, x, y, x_index, y_index):
+    length_2 = 0.08
+
+    pos_x_lower = x - length_2
+    pos_x_upper = x + length_2
+
+    pos_y_lower = y - length_2
+    pos_y_upper = y + length_2
+
+    x_lower = np.searchsorted(x_index, pos_x_lower)
+    x_upper = np.searchsorted(x_index, pos_x_upper)
+
+    y_lower = np.searchsorted(y_index, pos_y_lower)
+    y_upper = np.searchsorted(y_index, pos_y_upper)
+
+    return img[y_lower:y_upper, x_lower:x_upper]
+
+
 def main():
+    FFT_KWARGS = {"norm": matplotlib.colors.LogNorm(vmin=1), "cmap": "Greys"}
+    IMSHOW_WARGS = {"cmap": "Greys"}
+    SIZE = 601
+
     lat = Lattice(10, 10)
     maske = np.ones((10, 10), dtype=bool)
     x, y = lat.get_from_mask(maske)
     img = image_from_pos(x, y)
 
-    fig, [axs,axs2] = plt.subplots(2, 3)
+    img = padding(img, SIZE, SIZE)
-    axs[0].imshow(img)
+
+    fig, [axs, axs2] = plt.subplots(2, 3)
+    axs[0].imshow(img, **IMSHOW_WARGS)
     img = gaussian(img)
-    axs[1].imshow(img)
+    axs[1].imshow(img, **IMSHOW_WARGS)
-    plt.pause(.1)
+    plt.pause(0.1)
 
     freqx, freqy, intens = fft(img)
-    axs[2].imshow(intens, extent = (np.min(freqx), np.max(
+    axs[2].imshow(
-        freqx), np.min(freqy), np.max(freqy)))
+        intens,
-    
+        extent=(np.min(freqx), np.max(freqx), np.min(freqy), np.max(freqy)),
+        **FFT_KWARGS,
+    )
+
+    intens_rut = intens
     maske = np.zeros((10, 10), dtype=bool)
     x, y = lat.get_from_mask(maske)
     img = image_from_pos(x, y)
+    img = padding(img, SIZE, SIZE)
 
-    axs2[0].imshow(img)
+    axs2[0].imshow(img, **IMSHOW_WARGS)
     img = gaussian(img)
-    axs2[1].imshow(img)
+    axs2[1].imshow(img, **IMSHOW_WARGS)
-    plt.pause(.1)
+    plt.pause(0.1)
 
     freqx, freqy, intens = fft(img)
-    axs2[2].imshow(intens, extent = (np.min(freqx), np.max(
+    axs2[2].imshow(
-        freqx), np.min(freqy), np.max(freqy)))
+        intens,
+        extent=(np.min(freqx), np.max(freqx), np.min(freqy), np.max(freqy)),
+        **FFT_KWARGS,
+    )
+
+    # Create a Rectangle patch
+    # Add the patch to the Axes
+    point_x, point_y = reci_rutile()
+    for px, py in zip(point_x, point_y):
+        rect = rect_at_point(px, py, "r")
+        axs[2].add_patch(rect)
+        axs[2].text(
+            px, py, f"{np.sum(extract_rect(intens_rut, px, py, freqx, freqy)):0.2}"
+        )
+        rect = rect_at_point(px, py, "r")
+        axs2[2].add_patch(rect)
+        axs2[2].text(
+            px, py, f"{np.sum(extract_rect(intens, px, py, freqx, freqy)):0.2}"
+        )
+
+    point_x, point_y = reci_mono()
+    for px, py in zip(point_x, point_y):
+        # rect = rect_at_point(px, py,"b")
+        # axs[2].add_patch(rect)
+        rect = rect_at_point(px, py, "b")
+        axs2[2].add_patch(rect)
+        axs2[2].text(
+            px, py, f"{np.sum(extract_rect(intens, px, py, freqx, freqy)):0.2}"
+        )
+
+    axs[2].set_xlim(-1.0, 1.0)
+    axs[2].set_ylim(-1.0, 1.0)
+    axs2[2].set_xlim(-1.0, 1.0)
+    axs2[2].set_ylim(-1.0, 1.0)
+
+    plt.figure()
+    diff = intens_rut - intens
+    plt.imshow(
+        diff, extent=(np.min(freqx), np.max(freqx), np.min(freqy), np.max(freqy))
+    )
+
+    plt.xlim(-1.0, 1.0)
+    plt.ylim(-1.0, 1.0)
     plt.show()