added voronoi

2d_fourie/analysis.py

@@ -0,0 +1,81 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import glob
+
+
+def eval_data(file):
+    data = np.load(file)
+    percentage = data["percentage"]
+    out = data["out"]
+    out = np.array(out)
+    print(out.shape)
+    fig, all_axs = plt.subplots(2, out.shape[0])
+    axs = all_axs[0, :]
+    axs2 = all_axs[1, :]
+    for o, ax, ax2,  lab in zip(out, axs, axs2, ["rutile", "mono_twin", "mono"]):
+        # ax.plot(percentage, o/np.max(o, axis=0))
+        ax.plot(percentage, o/o[0])
+        # ax.plot(percentage, o)
+        o = np.mean(o, axis=1)
+        o = o/o[0]
+        ax2.plot(percentage, o)
+        ax2.plot([0, 1], [o[0], o[-1]], "k:")
+        ax.set_title(lab)
+    if "ising" in file:
+        fig.suptitle("Ising")
+    else:
+        fig.suptitle("Random")
+    fig.savefig(f"{file}.png")
+
+
+def parse_lists(out):
+    lists = []
+    for o in out:
+        lists.append(np.stack(o))
+
+    max = 0
+    for l in lists:
+        print(l.shape)
+        if max < l.shape[1]:
+            max = l.shape[1]
+    lists = [np.pad(l, ((0, 0), (0, max-l.shape[1]))) for l in lists]
+    for l in lists:
+        print(l.shape)
+    return np.stack(lists)
+
+
+def eval_data_print(file):
+    data = np.load(file, allow_pickle=True)
+    percentage = data["percentage"]
+    out = parse_lists(data["out"])
+    out = np.array(out)
+    print(out.shape)
+    out = out[[0, 2], :, :]
+    print(out.shape)
+    fig, all_axs = plt.subplots(2, out.shape[0])
+    axs = all_axs[0, :]
+    axs2 = all_axs[1, :]
+    for o, ax, ax2,  lab in zip(out, axs, axs2, ["rutile", "monoclinic", "mono"]):
+        # ax.plot(percentage, o/np.max(o, axis=0))
+        ax.plot(percentage, o/o[0])
+        # ax.plot(percentage, o)
+        o = np.mean(o, axis=1)
+        o = o/o[0]
+        ax2.plot(percentage, o)
+        ax2.plot([0, 1], [o[0], o[-1]], "k:")
+        ax.set_title(lab)
+    if "ising" in file:
+        fig.suptitle("Ising")
+    else:
+        fig.suptitle("Random")
+    for ax in all_axs.flatten():
+        ax.set_xlabel("Rutile Phase")
+        ax.set_ylabel("Normalized Intensity")
+
+    plt.tight_layout()
+
+
+if __name__ == "__main__":
+    for f in glob.glob("*.npz"):
+        eval_data_print(f)
+    plt.show()

2d_fourie/extractors.py

@@ -0,0 +1,195 @@
+import numpy as np
+from scipy.spatial import Voronoi
+import cv2
+
+
+class Image_Wrapper:
+    def __init__(self, img, x_lower, x_res, y_lower, y_res):
+        self.img = img
+        self.x_lower = x_lower
+        self.y_lower = y_lower
+        self.x_res = x_res
+        self.y_res = y_res
+
+        self.x_upper = self.x_lower + self.img.shape[0]*self.x_res - self.x_res
+        self.y_upper = self.y_lower + self.img.shape[1]*self.y_res - self.x_res
+
+    def __init__(self, img, fx, fy):
+        self.img = img
+        self.x_lower = np.min(fx)
+        self.y_lower = np.min(fy)
+        self.x_upper = np.max(fx)
+        self.y_upper = np.max(fy)
+
+        self.x_res = (self.x_upper - self.x_lower) / self.img.shape[0]
+        self.y_res = (self.y_upper - self.y_lower) / self.img.shape[0]
+
+    def val2pos(self, x, y):
+        x = (x - self.x_lower) / self.x_res
+        y = (y - self.y_lower) / self.y_res
+        return x, y
+
+    def check_bounds(self, xl, yl, xu, yu):
+        if xl > self.img.shape[0]:
+            print("xl lim")
+            return False
+        if yl > self.img.shape[1]:
+            print("yl lim")
+            return False
+        if xu < 0:
+            print("xu lim")
+            return False
+        if yu < 0:
+            print("yu lim")
+            return False
+        return True
+
+    def clean_bounds(self, xl, yl, xu, yu):
+        if xl < 0:
+            xl = 0
+        if yl < 0:
+            yl = 0
+
+        if xu > self.img.shape[0]:
+            xu = self.img.shape[0]
+        if yu > self.img.shape[1]:
+            yu = self.img.shape[1]
+
+        return xl, yl, xu, yu
+
+    def ext(self):
+        return [self.x_lower, self.x_upper, self.y_lower, self.y_upper]
+
+
+class Voronoi_Evaluator:
+    def __init__(self, points, eval_points):
+        self.eval_points = eval_points
+        self.vor = Voronoi(points)
+
+    def __init__(self, list_points):
+        points = np.concatenate(list_points, axis=0)
+        self.eval_points = []
+        start = 0
+        for l in list_points:
+            stop = l.shape[0]
+            self.eval_points.append(np.arange(start, start + stop))
+            start += stop
+        self.vor = Voronoi(points)
+
+    def extract(self, img: Image_Wrapper):
+        all = []
+        for ev_points in self.eval_points:
+            temp = []
+            region_mask = self.vor.point_region[ev_points]
+            for i in np.array(self.vor.regions)[region_mask]:
+                if -1 in i:
+                    print("Contains outside points")
+                    continue
+                if len(i) == 0:
+                    print("Contains outside points")
+                    continue
+                pts = self.vor.vertices[i]
+                pts = np.stack(img.val2pos(
+                    pts[:, 0], pts[:, 1])).astype(np.int32).T
+                mask = np.zeros_like(img.img)
+                cv2.fillConvexPoly(mask, pts, 1)
+                mask = mask > 0  # To convert to Boolean
+                temp.append(img.img[mask])
+                img.img[mask] = -1
+            all.append(temp)
+        return all
+
+    def extract_paint(self, img: Image_Wrapper):
+        counter  =  1
+        for ev_points in self.eval_points:
+            region_mask = self.vor.point_region[ev_points]
+            print(region_mask)
+            for i in np.array(self.vor.regions)[region_mask]:
+                if -1 in i:
+                    print("Contains outside points")
+                    continue
+                if len(i) == 0:
+                    print("Contains outside points")
+                    continue
+                pts = self.vor.vertices[i]
+                pts = np.stack(img.val2pos(
+                    pts[:, 0], pts[:, 1])).astype(np.int32).T
+                mask = np.zeros_like(img.img)
+                cv2.fillConvexPoly(mask, pts, 1)
+                mask = mask > 0  # To convert to Boolean
+                img.img[mask] = counter
+                counter += 1
+        return img.img
+
+
+class Rect_Evaluator:
+    def __init__(self, points, eval_points):
+        self.eval_points = eval_points
+        self.points = points
+        self.length = 4
+
+    def __init__(self, list_points):
+        self.points = np.concatenate(list_points, axis=0)
+        self.eval_points = []
+        start = 0
+        for l in list_points:
+            stop = l.shape[0]
+            self.eval_points.append(np.arange(start, start + stop))
+            start += stop
+        print(self.points.shape)
+        print(start, " from ")
+        self.length = 4
+
+    def extract(self, img: Image_Wrapper):
+        all = []
+        for ev_points in self.eval_points:
+            temp = []
+            for x, y in self.points[ev_points]:
+                x, y = img.val2pos(x, y)
+                x_lower = int(x - self.length)
+                y_lower = int(y - self.length)
+                x_upper = int(x + self.length + 1)
+                y_upper = int(y + self.length + 1)
+                if img.check_bounds(x_lower, y_lower, x_upper, y_upper):
+                    x_lower, y_lower, x_upper, y_upper = img.clean_bounds(
+                        x_lower, y_lower, x_upper, y_upper)
+                    temp.append(img.img[x_lower:x_upper])
+            all.append(temp)
+        return all
+
+    def extract_paint(self, img: Image_Wrapper):
+        val = np.nan
+        for ev_points in self.eval_points:
+            for x, y in self.points[ev_points]:
+                x, y = img.val2pos(x, y)
+                x_lower = int(x - self.length)
+                y_lower = int(y - self.length)
+                x_upper = int(x + self.length + 1)
+                y_upper = int(y + self.length + 1)
+                if img.check_bounds(x_lower, y_lower, x_upper, y_upper):
+                    x_lower, y_lower, x_upper, y_upper = img.clean_bounds(
+                        x_lower, y_lower, x_upper, y_upper)
+
+                    img.img[y_lower:y_upper, x_lower:x_upper] = val
+        return img.img
+
+#
+# def main():
+#     np.random.seed(10)
+#     points = (np.random.rand(100, 2)-0.5) * 2
+#     voro = Voronoi_Evaluator(points, [[1],[2]])
+#     rect = Rect_Evaluator(points, [[1], [2]])
+#     Z = np.ones((1000, 1000))
+#     img = Image_Wrapper(Z, -5, .01, -5, .01)
+#     voro.extract(img)
+#     rect.extract(img)
+#
+#     plt.scatter(points[[1], 0], points[[1], 1])
+#     plt.scatter(points[[2], 0], points[[2], 1])
+#     plt.imshow(img.img, extent=img.ext(), origin="lower")
+#     #plt.imshow(img.img, origin="lower")
+#     plt.show()
+#
+#
+# if __name__ == "__main__":
+#     main()

2d_fourie/lattices.py

@@ -18,6 +18,7 @@ class Lattice:
     def reci(self):
         pass
 
+
 class SCC_Lattice(Lattice):
     def __init__(self, x_len, y_len):
         x = np.arange(x_len) * 5
@@ -62,8 +63,6 @@ class VO2_Lattice(Lattice):
 
         offset_a_r, offset_c_r = self._mono_2_rutile(offset_c_m, offset_a_m)
 
-        print("A_r: ", offset_a_r, "C_r: ", offset_c_r)
-
         self.atom_x_mono = offset_a_r + X * \
             self.base_c_r + np.mod(Y, 4) * 0.5 * self.base_c_r
         self.atom_x_mono[np.mod(X, 2) == 0] -= 2 * offset_a_r
@@ -105,12 +104,11 @@ class VO2_Lattice(Lattice):
         return inplace_pos_x, inplace_pos_y
 
     def reci_rutile(self):
-       x = np.arange(-2, 3)
+        x = np.arange(-20, 21)
-       y = np.arange(-2, 3)
+        y = np.arange(-20, 21)
         X, Y = np.meshgrid(x, y)
         return (X * 0.22 + Y * 0.44).flatten(), (X * 0.349).flatten()
 
-
     def reci_mono(self):
         x, y = self.reci_rutile()
         return x + 0.1083, y + 0.1719
@@ -120,4 +118,14 @@ class VO2_Lattice(Lattice):
         return x - 0.1083, y + 0.1719
 
     def reci(self):
-        return [self.reci_rutile(), self.reci_mono(), self.reci_mono_2()]
+        cutoff = 5.
+        x, y = self.reci_rutile()
+        mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff)
+        p1 = (x[mask], y[mask])
+        x, y = self.reci_mono()
+        mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff)
+        p2 = (x[mask], y[mask])
+        x, y = self.reci_mono_2()
+        mask = np.logical_and(np.abs(x) < cutoff, np.abs(y) < cutoff)
+        p3 = (x[mask], y[mask])
+        return [p1, p2, p3]

2d_fourie/main.py

@@ -7,11 +7,13 @@ import matplotlib
 import scipy
 import scipy.signal
 import tqdm
+from extractors import Rect_Evaluator, Voronoi_Evaluator, Image_Wrapper
 
 
 class Plotter:
     def __init__(self, lat):
         self.lattice = lat
+        self.length_2 = 0.05
 
     def reduce(self, arr):
         arr = np.array(arr)
@@ -19,45 +21,8 @@ class Plotter:
         return np.mean(arr)
         # return np.sum(arr[np.argpartition(arr, -8)[-8:]])
 
-    def extract_rect(self, img, x, y, x_index, y_index):
-        length_2 = 0.01
-
-        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)
-
-        # fix different number of spins possible
-        if x_upper - x_lower < 10:
-            x_upper += 1
-        if y_upper - y_lower < 10:
-            y_upper += 1
-        return img[y_lower:y_upper, x_lower:x_upper]
-
-    def helper(self, ax, freqx, freqy, intens):
-        reci_lattice = self.lattice.reci()
-        for tup, col in zip(reci_lattice, ["r", "b", "g"]):
-            point_x, point_y = tup
-            point_x = point_x.flatten()
-            point_y = point_y.flatten()
-            for px, py in zip(point_x, point_y):
-                rect = self.rect_at_point(px, py, col)
-                ax.add_patch(rect)
-                sum = self.extract_rect(intens, px, py, freqx, freqy)
-                ax.text(
-                    px, py, f"{self.reduce(sum):2.2}", clip_on=True
-                )
-        return intens
-
     def rect_at_point(self, x, y, color):
-        length_2 = 0.01
+        length_2 = self.length_2
         rect = patches.Rectangle(
             (x - length_2, y - length_2),
             2 * length_2,
@@ -70,24 +35,21 @@ class Plotter:
 
     def plot(self, freqx, freqy, intens, ax_log=None, ax_lin=None, vmax=None):
         if ax_log:
-            intens = self.helper(ax_lin, freqx, freqy, intens)
             t = ax_log.imshow(
                 intens,
                 extent=(np.min(freqx), np.max(freqx),
                         np.min(freqy), np.max(freqy)),
-                norm=matplotlib.colors.LogNorm(vmin=10, vmax=vmax),
+                norm=matplotlib.colors.LogNorm(vmin=10),
                 cmap="viridis",
                 origin="lower"
             )
             plt.colorbar(t, ax=ax_log)
-            self.helper(ax_log, freqx, freqy, intens)
         if ax_lin:
-            intens = self.helper(ax_lin, freqx, freqy, intens)
             t = ax_lin.imshow(
                 intens,
                 extent=(np.min(freqx), np.max(freqx),
                         np.min(freqy), np.max(freqy)),
-                vmax=vmax,
+                #vmax=vmax,
                 cmap="viridis",
                 origin="lower"
             )
@@ -117,7 +79,6 @@ def test_square():
     plt.pause(0.1)
     fx, fy, intens = si.fft()
     plot.plot(fx, fy, intens, axs[1, 0], axs[1, 1])
-    print("Done")
     plt.savefig("test.png")
     plt.show()
 
@@ -126,6 +87,10 @@ def test_mixed():
     LEN = 40
     lat = VO2_Lattice(LEN, LEN)
     plot = Plotter(lat)
+    all_rutile = np.stack(lat.reci()[0]).T
+    all_mono = np.stack(lat.reci()[1]).T
+    all_mono2 = np.stack(lat.reci()[2]).T
+    rect = Voronoi_Evaluator([all_rutile, all_mono, all_mono2])
 
     pos_x, pos_y = lat.get_from_mask(np.zeros((40, 40)))
     si = SpinImage(pos_x, pos_y)
@@ -146,6 +111,13 @@ def test_mixed():
     si.pad_it_square(10)
     fx, fy, intens_mixed = si.fft()
 
+    intens_rutile = rect.extract_paint(
+        Image_Wrapper(intens_rutile, fx=fx, fy=fy))
+    intens_mono = rect.extract_paint(
+        Image_Wrapper(intens_mono, fx=fx, fy=fy))
+    intens_mixed = rect.extract_paint(
+        Image_Wrapper(intens_mixed, fx=fx, fy=fy))
+
     fig, axs = plt.subplots(2, 3)
     plot.plot(freqx=fx, freqy=fy, intens=intens_rutile,
               ax_log=axs[0, 0], ax_lin=axs[1, 0], vmax=10e7)
@@ -154,8 +126,6 @@ def test_mixed():
     plot.plot(freqx=fx, freqy=fy, intens=intens_mixed,
               ax_log=axs[0, 1], ax_lin=axs[1, 1], vmax=10e7)
 
-    print(np.sum(intens_mono), np.sum(intens_rutile), np.sum(intens_mixed))
-
     for ax in axs.flatten():
         ax.set_xlim(-1, 1)
         ax.set_ylim(-1, 1)
@@ -163,7 +133,8 @@ def test_mixed():
     plt.show()
 
 
-def random():
+def random(seed):
+    np.random.seed(seed)
     LEN = 40
     lat = VO2_Lattice(LEN, LEN)
     plot = Plotter(lat)
@@ -178,291 +149,97 @@ def random():
     for i in tqdm.tqdm(ind):
         maske[np.unravel_index(i, (LEN, LEN))] = True
         counter += 1
-        if np.mod(counter, 20) != 0:
+        if np.mod(counter, 100) != 0:
             continue
 
         pos_x, pos_y = lat.get_from_mask(maske)
         si = SpinImage(pos_x, pos_y)
         si.pad_it_square(10)
-        si.gaussian(LEN)
         fx, fy, intens = si.fft()
 
         for tup, lis in zip(reci_lattice, out):
             point_x, point_y = tup
             point_x = point_x.flatten()
             point_y = point_y.flatten()
-            sum = 0.
+            peaks = []
             for px, py in zip(point_x, point_y):
-                sum += np.sum(plot.extract_rect(intens, px, py, fx, fy))
+                peaks.append(np.sum(plot.extract_rect(intens, px, py, fx, fy)))
+            lis.append(peaks)
 
-            lis.append(sum)
+        percentage.append(np.sum(maske))
+
+    percentage = np.array(percentage)
+    percentage /= np.max(percentage)
+
+    np.savez(f"random_{seed}.npz", percentage=percentage, out=out)
+
+
+def sample_index(p):
+    i = np.random.choice(np.arange(p.size), p=p.ravel())
+    return np.unravel_index(i, p.shape)
+
+
+def ising(seed):
+    np.random.seed(seed)
+    LEN = 80
+    temp = 0.1
+    maske = np.zeros((LEN, LEN), dtype=bool)
+
+    lat = VO2_Lattice(LEN, LEN)
+    plot = Plotter(lat)
+
+    reci_lattice = lat.reci()
+    out = [[] for x in range(len(reci_lattice))]
+    percentage = []
+    counter = 0
+    for i in tqdm.tqdm(range(LEN*LEN)):
+        probability = np.roll(maske, 1, axis=0).astype(float)
+        probability += np.roll(maske, -1, axis=0).astype(float)
+        probability += np.roll(maske, 1, axis=1).astype(float)
+        probability += np.roll(maske, -1, axis=1).astype(float)
+
+        probability = np.exp(probability/temp)
+        probability[maske] = 0
+        probability /= np.sum(probability)
+        maske[sample_index(probability)] = True
+
+        counter += 1
+        if np.mod(counter, 100) != 0:
+            continue
+        pos_x, pos_y = lat.get_from_mask(maske)
+        si = SpinImage(pos_x, pos_y)
+        si.pad_it_square(10)
+        # si.gaussian(LEN)
+        fx, fy, intens = si.fft()
+
+        for tup, lis in zip(reci_lattice, out):
+            point_x, point_y = tup
+            point_x = point_x.flatten()
+            point_y = point_y.flatten()
+            peaks = []
+            for px, py in zip(point_x, point_y):
+                peaks.append(np.sum(plot.extract_rect(intens, px, py, fx, fy)))
+
+            lis.append(peaks)
         percentage.append(np.mean(maske))
+    percentage = np.array(percentage)
+    percentage /= np.max(percentage)
 
-    for o in out:
+    np.savez(f"ising_{temp}_{seed}.npz", percentage=percentage, out=out)
-        plt.scatter(percentage, o/o[0])
+    # for o in out:
-        plt.plot([0,1], [o[0], o[-1]])
+    #    plt.scatter(percentage, o/o[0])
-    plt.show()
+    #    plt.plot([0, 1], [1, o[-1]/o[0]])
+    # plt.pause(1)
 
 
 if __name__ == "__main__":
     # test_square()
-    # test_mixed()
+    test_mixed()
-    random()
+    plt.show()
-# def test_lattice():
+    # random()
-#    lat = VO2_Lattice(10, 10)
+    np.random.seed(1234)
-#    maske = np.zeros((10, 10), dtype=bool)
+    for i in np.random.randint(0, 10000, 2):
-#    x, y = lat.get_from_mask(maske)
+        random(i)
-#
+        ising(i)
-#    plt.scatter(x, y)
+
-#    maske = np.invert(maske)
+    plt.show()
-#    x, y = lat.get_from_mask(maske)
-#    plt.scatter(x, y)
-#
-#    maske[:3, :5] = False
-#    x, y = lat.get_from_mask(maske)
-#    plt.scatter(x, y)
-#    plt.show()
-#
-#
-# self.resolution = 0.1
-# CMAP = "Greys"
-#
-#
-# def test_img():
-#    lat = VO2_Lattice(10, 10)
-#    maske = np.ones((10, 10), dtype=bool)
-#    x, y = lat.get_from_mask(maske)
-#    img = image_from_pos(x, y)
-#    plt.imshow(img.T, origin="lower", extent=(0, np.max(x), 0, np.max(y)))
-#    plt.scatter(x, y)
-#    plt.show()
-#
-#
-# def gaussian(img):
-#    x = np.arange(-self.resolution * img.shape[0]/2,
-#                  self.resolution * img.shape[0]/2, self.resolution)
-#    y = np.arange(-self.resolution * img.shape[1]/2,
-#                  self.resolution * img.shape[1]/2, self.resolution)
-#    X, Y = np.meshgrid(x, y)
-#    sigma = self.resolution * img.shape[0] / 10
-#    print("Sigma: ", sigma)
-#    z = (
-#        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 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)
-#
-#    # fix different number of spins possible
-#    if x_upper - x_lower < 10:
-#        x_upper += 1
-#    if y_upper - y_lower < 10:
-#        y_upper += 1
-#
-#    return img[y_lower:y_upper, x_lower:x_upper]
-#
-#
-# def extract_peaks(freqx, freqy, intens):
-#    rutile = []
-#    point_x, point_y = reci_rutile()
-#    for px, py in zip(point_x, point_y):
-#        rutile.append(reduce(extract_rect(intens, px, py, freqx, freqy)))
-#
-#    mono = []
-#    point_x, point_y = reci_mono()
-#    for px, py in zip(point_x, point_y):
-#        mono.append(reduce(extract_rect(intens, px, py, freqx, freqy)))
-#    return rutile, mono
-#
-#
-# def plot(ax, freqx, freqy, intens):
-#    point_x, point_y = reci_rutile()
-#    for px, py in zip(point_x, point_y):
-#        rect = rect_at_point(px, py, "r")
-#        ax.add_patch(rect)
-#        ax.text(
-#            px, py, f"{reduce(extract_rect(intens, px, py, freqx, freqy)):2.2}", clip_on=True
-#        )
-#
-#    point_x, point_y = reci_mono()
-#    for px, py in zip(point_x, point_y):
-#        rect = rect_at_point(px, py, "b")
-#        ax.add_patch(rect)
-#        ax.text(
-#            px, py, f"{reduce(extract_rect(intens, px, py, freqx, freqy)):2.2}", clip_on=True
-#        )
-#    ax.imshow(
-#        intens,
-#        extent=(np.min(freqx), np.max(freqx), np.min(freqy), np.max(freqy)),
-#        norm=matplotlib.colors.LogNorm(),
-#        cmap="Greys"
-#    )
-#
-#
-# def test_all():
-#    LEN = 100
-#    SIZE = 60 * LEN + 1
-#    quad = np.ones((3, 3))
-#
-#    fig, ax = plt.subplots(1, 3)
-#
-#    lat = VO2_Lattice(LEN, LEN)
-#    maske = np.ones((LEN, LEN), dtype=bool)
-#    x, y = lat.get_from_mask(maske)
-#    img = image_from_pos(x, y)
-#    img = padding(img, SIZE, SIZE)
-#    #img = scipy.signal.convolve2d(img, quad)
-#    img = gaussian(img)
-#    freqx, freqy, intens_rutile = fft(img)
-#
-#    img = scipy.signal.convolve2d(img, quad)
-#    ax[0].imshow(img)
-#
-#    maske = np.zeros((LEN, LEN), dtype=bool)
-#    x, y = lat.get_from_mask(maske)
-#    img = image_from_pos(x, y)
-#    img = padding(img, SIZE, SIZE)
-#    img = gaussian(img)
-#    freqx, freqy, intens_mono = fft(img)
-#
-#    img = scipy.signal.convolve2d(img, quad)
-#    ax[2].imshow(img)
-#
-#    maske = np.zeros((LEN, LEN), dtype=bool)
-#    ind = np.arange(LEN*LEN)
-#    np.random.shuffle(ind)
-#    ind = np.unravel_index(ind[:int(LEN*LEN/2)], (LEN, LEN))
-#    maske[ind] = True
-#    x, y = lat.get_from_mask(maske)
-#    img = image_from_pos(x, y)
-#    img = padding(img, SIZE, SIZE)
-#    img = gaussian(img)
-#    freqx, freqy, intens_mono = fft(img)
-#
-#    img = scipy.signal.convolve2d(img, quad)
-#    ax[1].imshow(img)
-#
-#    print(np.mean(maske))
-#    x, y = lat.get_from_mask(maske)
-#    img = image_from_pos(x, y)
-#    img = padding(img, SIZE, SIZE)
-#    img = gaussian(img)
-#    freqx, freqy, intens_50 = fft(img)
-#
-#    fig, axs = plt.subplots(1, 3)
-#    plot(axs[0], freqx=freqx, freqy=freqy, intens=intens_rutile)
-#    plot(axs[2], freqx=freqx, freqy=freqy, intens=intens_mono)
-#    plot(axs[1], freqx=freqx, freqy=freqy, intens=intens_50)
-#    axs[0].set_title("Rutile")
-#    axs[2].set_title("Mono")
-#    axs[1].set_title("50/50")
-#
-#    for ax in axs:
-#        ax.set_xlim(-1.0, 1.0)
-#        ax.set_ylim(-1.0, 1.0)
-#
-#
-# def eval(maske, lat, LEN):
-#    x, y = lat.get_from_mask(maske)
-#    SIZE = 60 * LEN + 1
-#    img = image_from_pos(x, y)
-#    img = padding(img, SIZE, SIZE)
-#    img = gaussian(img)
-#    freqx, freqy, intens = fft(img)
-#    return extract_peaks(freqx, freqy, intens)
-#
-#
-# def reduce(arr):
-#    arr = np.array(arr)
-#    arr = arr.flatten()
-#    return np.sum(arr[np.argpartition(arr, -8)[-8:]])
-#
-#
-# def main():
-#    LEN = 80
-#    lat = VO2_Lattice(LEN, LEN)
-#    maske = np.zeros((LEN, LEN), dtype=bool)
-#    ind = np.arange(LEN*LEN)
-#    np.random.shuffle(ind)
-#    percentage = []
-#    rutile = []
-#    monoclinic = []
-#    counter = 0
-#    for i in tqdm.tqdm(ind):
-#        i_unravel = np.unravel_index(i, (LEN, LEN))
-#        maske[i_unravel] = True
-#        if np.mod(counter, 300) == 0:
-#            rut, mono = eval(maske, lat, LEN)
-#            percentage.append(np.mean(maske))
-#            rutile.append(reduce(rut))
-#            monoclinic.append(reduce(mono))
-#        counter += 1
-#
-#    print(len(percentage), len(mono), len(rutile))
-#    print(mono)
-#    plt.figure()
-#    plt.scatter(percentage, np.array(monoclinic)/monoclinic[0], label="mono")
-#    plt.scatter(percentage, np.array(rutile)/rutile[0], label="rut")
-#    plt.legend()
-#
-#
-# if __name__ == "__main__":
-#    test_all()
-#    # main()
-#    plt.show()

2d_fourie/plot.py

@@ -0,0 +1,5 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+if __name__ == "__main__":
+    pass

2d_fourie/spin_image.py

@@ -1,6 +1,6 @@
 import numpy as np
 import scipy.fftpack as sfft
-
+import scipy
 
 class SpinImage:
     resolution = 0.1
@@ -33,12 +33,10 @@ class SpinImage:
     def pad_it_square(self, additional_pad=0):
         h = self.img.shape[0]
         w = self.img.shape[1]
-        print(h, w)
         xx = np.maximum(h, w) + 2 * additional_pad
         yy = xx
         self.length_x = xx * self.resolution
         self.length_y = yy * self.resolution
-        print("Pad to: ", xx, yy)
 
         a = (xx - h) // 2
         aa = xx - a - h
@@ -49,6 +47,18 @@ class SpinImage:
         self.img = np.pad(self.img, pad_width=(
             (a, aa), (b, bb)), mode="constant")
 
+    def percentage_gaussian(self, mask, sigma):
+        x = np.linspace(-self.length_x / 2,
+                      self.length_x / 2, mask.shape[0])
+        y = np.linspace(-self.length_y / 2,
+                      self.length_y / 2, mask.shape[1])
+        X, Y = np.meshgrid(x, y)
+        z = (
+            1 / (2 * np.pi * sigma * sigma)
+            * np.exp(-(X ** 2 / (2 * sigma ** 2) + Y ** 2 / (2 * sigma ** 2)))
+        )
+        return np.multiply(mask, z.T)
+
     def gaussian(self, sigma):
         x = np.arange(-self.length_x / 2,
                       self.length_x / 2, self.resolution)

2d_fourie/test_voronoi.py

@@ -0,0 +1,46 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.spatial import Voronoi, voronoi_plot_2d
+import cv2
+
+
+
+class Voronoi_Evaulator:
+    def __init__(self, points, eval_points):
+        self.eval_points = eval_points
+        self.vor = Voronoi(points)
+
+    def extract(self, Z):
+        for debug_num, ev_points in zip([-10, -5], self.eval_points):
+            region_mask = self.vor.point_region[ev_points]
+            print(region_mask)
+            for i in np.array(self.vor.regions)[region_mask]:
+                if -1 in i:
+                    print("Containse outside points")
+                    continue
+                if len(i) == 0:
+                    print("Containse outside points")
+                    continue
+                print(i)
+                pts = self.vor.vertices[i]
+                pts = (pts * 100).astype(np.int32)
+                print(pts)
+                mask = np.zeros((Z.shape[0], Z.shape[1]))
+                cv2.fillConvexPoly(mask, pts, 1)
+                mask = mask > 0  # To convert to Boolean
+                Z[mask] = debug_num
+        return Z
+
+if __name__ == "__main__":
+    np.random.seed(20)
+    points = (np.random.rand(100, 2)-0.1) * 2
+    voro = Voronoi_Evaulator(points, [[1, 4, 5], [2, 3, 6]])
+
+    x = np.linspace(0, 1, 100)
+    y = np.linspace(0, 1, 200)
+    X, Y = np.meshgrid(x, y)
+    Z = X*2 + Y
+
+    Z = voro.extract(Z)
+    plt.imshow(Z)
+    plt.show()