diff --git a/2d_fourie/extractors.py b/2d_fourie/extractors.py index f296d71..5b33452 100644 --- a/2d_fourie/extractors.py +++ b/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): - self.img = img - self.x_lower = x_lower - self.y_lower = y_lower - self.x_res = x_res - self.y_res = y_res + # 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 + # 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 @@ -61,11 +61,41 @@ class Image_Wrapper: 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) +class Evaluator: + def extract(self, img: Image_Wrapper): + all_val = [] + 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): - 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 generate_mask(self, img: Image_Wrapper): + self.mask = np.full_like(img.img, -1) - 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 + counter = -1 + region_mask = self.vor.point_region + for i in np.array(self.vor.regions, dtype=list)[region_mask]: + counter += 1 + if -1 in i: + continue + if len(i) == 0: + continue + pts = self.vor.vertices[i] + pts = np.stack(img.val2pos( + pts[:, 0], pts[:, 1])).astype(np.int32).T + if np.any(pts < 0): + continue + mask = np.zeros_like(img.img) + cv2.fillConvexPoly(mask, pts, 1) + mask = mask > 0 # To convert to Boolean + self.mask[mask] = counter -class Rect_Evaluator: - def __init__(self, points, eval_points): - self.eval_points = eval_points - self.points = points - self.length = 4 +class Rect_Evaluator(Evaluator): + # def __init__(self, points, eval_points): + # self.eval_points = eval_points + # self.points = points + # 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 = 4 + self.length = length - 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) + def generate_mask(self, img: Image_Wrapper): + self.mask = np.full_like(img.img, -1) + count = 0 + for x, y in self.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) + self.mask[y_lower:y_upper, x_lower:x_upper] = count + count += 1 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) diff --git a/2d_fourie/lattices.py b/2d_fourie/lattices.py index b3b17cf..cc139ce 100644 --- a/2d_fourie/lattices.py +++ b/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 - y = np.arange(-3, 3) * 0.2 + x = np.arange(-3, 4) * 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) - y = np.arange(-20, 21) + num = 20 + 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() diff --git a/2d_fourie/main.py b/2d_fourie/main.py index 2c18e14..cc3065e 100644 --- a/2d_fourie/main.py +++ b/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) - fx, fy, intens_rutile = si.fft() + si.pad_it_square(10, size=2300) + 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) - fx, fy, intens_mono = si.fft() + si.pad_it_square(10, size=2300) + 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( - 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(1, 2) + axs[0].imshow(rect.get_mask(), extent=img.ext(), origin="lower") + axs[0].plot(all_rutile[:, 0], all_rutile[:, 1], ".") + axs[1].plot(all_rutile[:, 0], all_rutile[:, 1], ".") + axs[0].plot(all_mono[:, 0], all_mono[:, 1], ".") + 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): - 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) - + iv, vv = voro.extract(img) + ir, vr = rect.extract(img) + for lis, val in zip(out_rect, vr): + lis.append(val) + for lis, val in zip(out_voro, vv): + lis.append(val) 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 = [[] 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(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.gaussian(LEN) + 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): - 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))) + iv, vv = voro.extract(img) + ir, vr = rect.extract(img) + for lis, val in zip(out_rect, vr): + lis.append(val) + for lis, val in zip(out_voro, vv): + lis.append(val) + percentage.append(np.sum(maske)) - lis.append(peaks) - percentage.append(np.mean(maske)) - percentage = np.array(percentage) + percentage = np.array(percentage, dtype=np.float64) percentage /= np.max(percentage) - np.savez(f"ising_{temp}_{seed}.npz", percentage=percentage, out=out) - # for o in out: - # plt.scatter(percentage, o/o[0]) - # plt.plot([0, 1], [1, o[-1]/o[0]]) - # plt.pause(1) + np.savez(f"ising_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"ising_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]) 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) diff --git a/2d_fourie/spin_image.py b/2d_fourie/spin_image.py index 6980f08..e7ecbef 100644 --- a/2d_fourie/spin_image.py +++ b/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)