Added voronoi support

2d_fourie/extractors.py

@@ -4,15 +4,15 @@ import cv2
 
 
 class Image_Wrapper:
-    def __init__(self, img, x_lower, x_res, y_lower, y_res):
+    # def __init__(self, img, x_lower, x_res, y_lower, y_res):
-        self.img = img
+    #    self.img = img
-        self.x_lower = x_lower
+    #    self.x_lower = x_lower
-        self.y_lower = y_lower
+    #    self.y_lower = y_lower
-        self.x_res = x_res
+    #    self.x_res = x_res
-        self.y_res = y_res
+    #    self.y_res = y_res
 
-        self.x_upper = self.x_lower + self.img.shape[0]*self.x_res - self.x_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
+    #    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
@@ -61,11 +61,41 @@ class Image_Wrapper:
         return [self.x_lower, self.x_upper, self.y_lower, self.y_upper]
 
 
-class Voronoi_Evaluator:
+class Evaluator:
-    def __init__(self, points, eval_points):
+    def extract(self, img: Image_Wrapper):
-        self.eval_points = eval_points
+        all_val = []
-        self.vor = Voronoi(points)
+        all_idx = []
+        for ev_points in self.eval_points:
+            temp_val = []
+            temp_idx = []
+            for num in ev_points:
+                if np.sum(self.mask == num) == 0:
+                    # print("Empty: ",num)
+                    continue
+                temp_val.append(np.mean(img.img[self.mask == num]))
+                temp_idx.append(num)
+            all_val.append(temp_val)
+            all_idx.append(temp_idx)
+        all_val.append([np.mean(img.img[self.mask == -1])])
+        all_idx.append([-1])
+        return all_idx, all_val
 
+    def debug(self, img: Image_Wrapper):
+        count = 1
+        for ev_points in self.eval_points:
+            if count != 1 and False:
+                count += 1
+                continue
+            for num in ev_points:
+                img.img[self.mask == num] += 10000 * num
+            count += 1
+        return img.img
+
+    def get_mask(self):
+        return self.mask
+
+
+class Voronoi_Evaluator(Evaluator):
     def __init__(self, list_points):
         points = np.concatenate(list_points, axis=0)
         self.eval_points = []
@@ -76,59 +106,35 @@ class Voronoi_Evaluator:
             start += stop
         self.vor = Voronoi(points)
 
-    def extract(self, img: Image_Wrapper):
+    def generate_mask(self, img: Image_Wrapper):
-        all = []
+        self.mask = np.full_like(img.img, -1)
-        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
-        counter  =  1
+        region_mask = self.vor.point_region
-        for ev_points in self.eval_points:
+        for i in np.array(self.vor.regions, dtype=list)[region_mask]:
-            region_mask = self.vor.point_region[ev_points]
+            counter += 1
-            print(region_mask)
+            if -1 in i:
-            for i in np.array(self.vor.regions)[region_mask]:
+                continue
-                if -1 in i:
+            if len(i) == 0:
-                    print("Contains outside points")
+                continue
-                    continue
+            pts = self.vor.vertices[i]
-                if len(i) == 0:
+            pts = np.stack(img.val2pos(
-                    print("Contains outside points")
+                pts[:, 0], pts[:, 1])).astype(np.int32).T
-                    continue
+            if np.any(pts < 0):
-                pts = self.vor.vertices[i]
+                continue
-                pts = np.stack(img.val2pos(
+            mask = np.zeros_like(img.img)
-                    pts[:, 0], pts[:, 1])).astype(np.int32).T
+            cv2.fillConvexPoly(mask, pts, 1)
-                mask = np.zeros_like(img.img)
+            mask = mask > 0  # To convert to Boolean
-                cv2.fillConvexPoly(mask, pts, 1)
+            self.mask[mask] = counter
-                mask = mask > 0  # To convert to Boolean
-                img.img[mask] = counter
-                counter += 1
-        return img.img
 
 
-class Rect_Evaluator:
+class Rect_Evaluator(Evaluator):
-    def __init__(self, points, eval_points):
+    # def __init__(self, points, eval_points):
-        self.eval_points = eval_points
+    #    self.eval_points = eval_points
-        self.points = points
+    #    self.points = points
-        self.length = 4
+    #    self.length = 4
 
-    def __init__(self, list_points):
+    def __init__(self, list_points, length=4):
         self.points = np.concatenate(list_points, axis=0)
         self.eval_points = []
         start = 0
@@ -136,43 +142,26 @@ class Rect_Evaluator:
             stop = l.shape[0]
             self.eval_points.append(np.arange(start, start + stop))
             start += stop
+            print(start)
         print(self.points.shape)
-        print(start, " from ")
+        self.length = length
-        self.length = 4
 
-    def extract(self, img: Image_Wrapper):
+    def generate_mask(self, img: Image_Wrapper):
-        all = []
+        self.mask = np.full_like(img.img, -1)
-        for ev_points in self.eval_points:
+        count = 0
-            temp = []
+        for x, y in self.points:
-            for x, y in self.points[ev_points]:
+            x, y = img.val2pos(x, y)
-                x, y = img.val2pos(x, y)
+            x_lower = int(x - self.length)
-                x_lower = int(x - self.length)
+            y_lower = int(y - self.length)
-                y_lower = int(y - self.length)
+            x_upper = int(x + self.length + 1)
-                x_upper = int(x + self.length + 1)
+            y_upper = int(y + self.length + 1)
-                y_upper = int(y + self.length + 1)
+            if img.check_bounds(x_lower, y_lower, x_upper, y_upper):
-                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.clean_bounds(
+                    x_lower, y_lower, x_upper, y_upper)
-                        x_lower, y_lower, x_upper, y_upper)
+                self.mask[y_lower:y_upper, x_lower:x_upper] = count
-                    temp.append(img.img[x_lower:x_upper])
+                count += 1
-            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)

2d_fourie/lattices.py

@@ -29,8 +29,8 @@ class SCC_Lattice(Lattice):
         return self.X, self.Y
 
     def reci(self):
-        x = np.arange(-3, 3) * 0.2
+        x = np.arange(-3, 4) * 0.2
-        y = np.arange(-3, 3) * 0.2
+        y = np.arange(-3, 4) * 0.2
         X, Y = np.meshgrid(x, y)
         return [(X, Y)]
 
@@ -104,8 +104,10 @@ class VO2_Lattice(Lattice):
         return inplace_pos_x, inplace_pos_y
 
     def reci_rutile(self):
-        x = np.arange(-20, 21)
+        num = 20
-        y = np.arange(-20, 21)
+        num = 2
+        x = np.arange(-num, num + 1)
+        y = np.arange(-num, num + 1)
         X, Y = np.meshgrid(x, y)
         return (X * 0.22 + Y * 0.44).flatten(), (X * 0.349).flatten()
 

2d_fourie/main.py

@@ -49,7 +49,7 @@ class Plotter:
                 intens,
                 extent=(np.min(freqx), np.max(freqx),
                         np.min(freqy), np.max(freqy)),
-                #vmax=vmax,
+                # vmax=vmax,
                 cmap="viridis",
                 origin="lower"
             )
@@ -83,6 +83,20 @@ def test_square():
     plt.show()
 
 
+def helper(intens, fx, fy, voro, rect):
+    print(np.mean(intens))
+    v_idx, v_val = voro.extract(Image_Wrapper(intens, fx, fy))
+    r_idx, r_val = rect.extract(Image_Wrapper(intens, fx, fy))
+
+    for iv, av, ir, ar in zip(v_idx, v_val, r_idx, r_val):
+        av = np.array(av)
+        ar = np.array(ar)
+        mask = np.isin(ir, iv)
+        print("Test: ", np.all(np.isclose(
+            ar[mask], av)), np.isclose(ar[mask], av))
+        print(iv, "\n",  ar[mask], "\n", av, "\n\n\n")
+
+
 def test_mixed():
     LEN = 40
     lat = VO2_Lattice(LEN, LEN)
@@ -90,17 +104,20 @@ def test_mixed():
     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])
+    voro = Voronoi_Evaluator([all_rutile, all_mono, all_mono2])
+    rect = Rect_Evaluator([all_rutile, all_mono, all_mono2])
+    #voro = Voronoi_Evaluator([all_mono])
+    #rect = Rect_Evaluator([all_mono])
 
     pos_x, pos_y = lat.get_from_mask(np.zeros((40, 40)))
     si = SpinImage(pos_x, pos_y)
-    si.pad_it_square(10)
+    si.pad_it_square(10, size=2300)
-    fx, fy, intens_rutile = si.fft()
+    fx, fy, intens_mono = si.fft()
 
     pos_x, pos_y = lat.get_from_mask(np.ones((40, 40)))
     si = SpinImage(pos_x, pos_y)
-    si.pad_it_square(10)
+    si.pad_it_square(10, size=2300)
-    fx, fy, intens_mono = si.fft()
+    fx, fy, intens_rutile = si.fft()
 
     mask_misk = np.ones((40, 40))
     ind = np.arange(mask_misk.size)
@@ -108,44 +125,80 @@ def test_mixed():
     mask_misk[np.unravel_index(ind[:800], (40, 40))] = False
     pos_x, pos_y = lat.get_from_mask(mask_misk)
     si = SpinImage(pos_x, pos_y)
-    si.pad_it_square(10)
+    si.pad_it_square(10, size=2300)
     fx, fy, intens_mixed = si.fft()
+    img = Image_Wrapper(intens_mono, fx, fy)
+    voro.generate_mask(img)
+    rect.generate_mask(Image_Wrapper(intens_mono, fx, fy))
 
-    intens_rutile = rect.extract_paint(
+    fig, axs = plt.subplots(1, 2)
-        Image_Wrapper(intens_rutile, fx=fx, fy=fy))
+    axs[0].imshow(rect.get_mask(), extent=img.ext(), origin="lower")
-    intens_mono = rect.extract_paint(
+    axs[0].plot(all_rutile[:, 0], all_rutile[:, 1], ".")
-        Image_Wrapper(intens_mono, fx=fx, fy=fy))
+    axs[1].plot(all_rutile[:, 0], all_rutile[:, 1], ".")
-    intens_mixed = rect.extract_paint(
+    axs[0].plot(all_mono[:, 0], all_mono[:, 1], ".")
-        Image_Wrapper(intens_mixed, fx=fx, fy=fy))
+    axs[1].plot(all_mono[:, 0], all_mono[:, 1], ".")
+    axs[0].plot(all_mono2[:, 0], all_mono2[:, 1], ".")
+    axs[1].plot(all_mono2[:, 0], all_mono2[:, 1], ".")
+    axs[1].imshow(voro.get_mask(), extent=img.ext(), origin="lower")
+
+    print("mono")
+    helper(intens=intens_mono, fx=fx, fy=fy, voro=voro, rect=rect)
+    print("mixed")
+    helper(intens=intens_mixed, fx=fx, fy=fy, voro=voro, rect=rect)
+    print("rutile")
+    helper(intens=intens_rutile, fx=fx, fy=fy, voro=voro, rect=rect)
+
+    new_intens_mono = rect.debug(Image_Wrapper(intens_mono, fx, fy))
+    new_intens_mixed = rect.debug(Image_Wrapper(intens_mixed, fx, fy))
+    new_intens_rutile = rect.debug(Image_Wrapper(intens_rutile, fx, fy))
 
     fig, axs = plt.subplots(2, 3)
-    plot.plot(freqx=fx, freqy=fy, intens=intens_rutile,
+    plot.plot(freqx=fx, freqy=fy, intens=new_intens_rutile,
               ax_log=axs[0, 0], ax_lin=axs[1, 0], vmax=10e7)
-    plot.plot(freqx=fx, freqy=fy, intens=intens_mono,
+    plot.plot(freqx=fx, freqy=fy, intens=new_intens_mono,
               ax_log=axs[0, 2], ax_lin=axs[1, 2], vmax=10e7)
-    plot.plot(freqx=fx, freqy=fy, intens=intens_mixed,
+    plot.plot(freqx=fx, freqy=fy, intens=new_intens_mixed,
               ax_log=axs[0, 1], ax_lin=axs[1, 1], vmax=10e7)
 
     for ax in axs.flatten():
         ax.set_xlim(-1, 1)
         ax.set_ylim(-1, 1)
 
-    plt.show()
+    new_intens_mono = voro.debug(Image_Wrapper(new_intens_mono, fx, fy))
+    new_intens_mixed = voro.debug(Image_Wrapper(intens_mixed, fx, fy))
+    new_intens_rutile = voro.debug(Image_Wrapper(intens_rutile, fx, fy))
+
+    fig, axs = plt.subplots(2, 3)
+    plot.plot(freqx=fx, freqy=fy, intens=new_intens_rutile,
+              ax_log=axs[0, 0], ax_lin=axs[1, 0], vmax=10e7)
+    plot.plot(freqx=fx, freqy=fy, intens=new_intens_mono,
+              ax_log=axs[0, 2], ax_lin=axs[1, 2], vmax=10e7)
+    plot.plot(freqx=fx, freqy=fy, intens=new_intens_mixed,
+              ax_log=axs[0, 1], ax_lin=axs[1, 1], vmax=10e7)
+
+    for ax in axs.flatten():
+        ax.set_xlim(-1, 1)
+        ax.set_ylim(-1, 1)
 
 
 def random(seed):
     np.random.seed(seed)
     LEN = 40
     lat = VO2_Lattice(LEN, LEN)
-    plot = Plotter(lat)
     maske = np.zeros((LEN, LEN))
     ind = np.arange(LEN * LEN)
     np.random.shuffle(ind)
-    reci_lattice = lat.reci()
+    all_rutile = np.stack(lat.reci()[0]).T
+    all_mono = np.stack(lat.reci()[1]).T
+    all_mono2 = np.stack(lat.reci()[2]).T
+    voro = Voronoi_Evaluator([all_rutile, all_mono, all_mono2])
+    rect = Rect_Evaluator([all_rutile, all_mono, all_mono2])
 
-    out = [[] for x in range(len(reci_lattice))]
+    out_rect = [[] for x in range(len(lat.reci())+1)]
+    out_voro = [[] for x in range(len(lat.reci())+1)]
     percentage = []
     counter = 0
+    already_inited = False
     for i in tqdm.tqdm(ind):
         maske[np.unravel_index(i, (LEN, LEN))] = True
         counter += 1
@@ -154,24 +207,31 @@ def random(seed):
 
         pos_x, pos_y = lat.get_from_mask(maske)
         si = SpinImage(pos_x, pos_y)
-        si.pad_it_square(10)
+        si.pad_it_square(10, size=2300)
         fx, fy, intens = si.fft()
+        img = Image_Wrapper(intens, fx, fy)
+        if not already_inited:
+            voro.generate_mask(img)
+            rect.generate_mask(img)
+            already_inited = True
 
-        for tup, lis in zip(reci_lattice, out):
+        iv, vv = voro.extract(img)
-            point_x, point_y = tup
+        ir, vr = rect.extract(img)
-            point_x = point_x.flatten()
+        for lis, val in zip(out_rect, vr):
-            point_y = point_y.flatten()
+            lis.append(val)
-            peaks = []
+        for lis, val in zip(out_voro, vv):
-            for px, py in zip(point_x, point_y):
+            lis.append(val)
-                peaks.append(np.sum(plot.extract_rect(intens, px, py, fx, fy)))
-            lis.append(peaks)
-
         percentage.append(np.sum(maske))
 
     percentage = np.array(percentage)
     percentage /= np.max(percentage)
 
-    np.savez(f"random_{seed}.npz", percentage=percentage, out=out)
+    np.savez(f"random_rect_{seed}.npz",
+             percentage=percentage, out_1=out_rect[0],
+             out_2=out_rect[1], out_3=out_rect[2], out_4=out_rect[3])
+    np.savez(f"random_voro_{seed}.npz",
+             percentage=percentage, out_1=out_voro[0],
+             out_2=out_voro[1], out_3=out_voro[2], out_4=out_voro[3])
 
 
 def sample_index(p):
@@ -181,17 +241,22 @@ def sample_index(p):
 
 def ising(seed):
     np.random.seed(seed)
-    LEN = 80
+    LEN = 40
     temp = 0.1
     maske = np.zeros((LEN, LEN), dtype=bool)
 
     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
+    voro = Voronoi_Evaluator([all_rutile, all_mono, all_mono2])
+    rect = Rect_Evaluator([all_rutile, all_mono, all_mono2])
 
-    reci_lattice = lat.reci()
+    out_rect = [[] for x in range(len(lat.reci())+1)]
-    out = [[] for x in range(len(reci_lattice))]
+    out_voro = [[] for x in range(len(lat.reci())+1)]
     percentage = []
     counter = 0
+    already_inited = False
     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)
@@ -206,30 +271,34 @@ def ising(seed):
         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.pad_it_square(10, size=2300)
-        # si.gaussian(LEN)
         fx, fy, intens = si.fft()
+        img = Image_Wrapper(intens, fx, fy)
+        if not already_inited:
+            voro.generate_mask(img)
+            rect.generate_mask(img)
+            already_inited = True
 
-        for tup, lis in zip(reci_lattice, out):
+        iv, vv = voro.extract(img)
-            point_x, point_y = tup
+        ir, vr = rect.extract(img)
-            point_x = point_x.flatten()
+        for lis, val in zip(out_rect, vr):
-            point_y = point_y.flatten()
+            lis.append(val)
-            peaks = []
+        for lis, val in zip(out_voro, vv):
-            for px, py in zip(point_x, point_y):
+            lis.append(val)
-                peaks.append(np.sum(plot.extract_rect(intens, px, py, fx, fy)))
+        percentage.append(np.sum(maske))
 
-            lis.append(peaks)
+    percentage = np.array(percentage, dtype=np.float64)
-        percentage.append(np.mean(maske))
-    percentage = np.array(percentage)
     percentage /= np.max(percentage)
 
-    np.savez(f"ising_{temp}_{seed}.npz", percentage=percentage, out=out)
+    np.savez(f"ising_rect_{seed}.npz",
-    # for o in out:
+             percentage=percentage, out_1=out_rect[0],
-    #    plt.scatter(percentage, o/o[0])
+             out_2=out_rect[1], out_3=out_rect[2], out_4=out_rect[3])
-    #    plt.plot([0, 1], [1, o[-1]/o[0]])
+    np.savez(f"ising_voro_{seed}.npz",
-    # plt.pause(1)
+             percentage=percentage, out_1=out_voro[0],
+             out_2=out_voro[1], out_3=out_voro[2], out_4=out_voro[3])
 
 
 if __name__ == "__main__":
@@ -238,7 +307,7 @@ if __name__ == "__main__":
     plt.show()
     # random()
     np.random.seed(1234)
-    for i in np.random.randint(0, 10000, 2):
+    for i in np.random.randint(0, 10000, 1):
         random(i)
         ising(i)
 

2d_fourie/spin_image.py

@@ -2,6 +2,7 @@ import numpy as np
 import scipy.fftpack as sfft
 import scipy
 
+
 class SpinImage:
     resolution = 0.1
 
@@ -30,10 +31,12 @@ class SpinImage:
         fft_freqy_clean = sfft.fftshift(fft_freqy)
         return fft_freqx_clean, fft_freqy_clean, np.abs(Z_shift) ** 2
 
-    def pad_it_square(self, additional_pad=0):
+    def pad_it_square(self, additional_pad=0, size=None):
         h = self.img.shape[0]
         w = self.img.shape[1]
         xx = np.maximum(h, w) + 2 * additional_pad
+        if size is not None:
+            xx = np.maximum(xx, size)
         yy = xx
         self.length_x = xx * self.resolution
         self.length_y = yy * self.resolution
@@ -49,9 +52,9 @@ class SpinImage:
 
     def percentage_gaussian(self, mask, sigma):
         x = np.linspace(-self.length_x / 2,
-                      self.length_x / 2, mask.shape[0])
+                        self.length_x / 2, mask.shape[0])
         y = np.linspace(-self.length_y / 2,
-                      self.length_y / 2, mask.shape[1])
+                        self.length_y / 2, mask.shape[1])
         X, Y = np.meshgrid(x, y)
         z = (
             1 / (2 * np.pi * sigma * sigma)