From 4d379426499b3ab6394da38e98c3d9409813ec5f Mon Sep 17 00:00:00 2001 From: Jacob Date: Wed, 1 Mar 2023 12:50:53 +0100 Subject: [PATCH] added voronoi --- 2d_fourie/analysis.py | 81 ++++++++ 2d_fourie/extractors.py | 195 ++++++++++++++++++ 2d_fourie/lattices.py | 40 ++-- 2d_fourie/main.py | 413 +++++++++----------------------------- 2d_fourie/plot.py | 5 + 2d_fourie/spin_image.py | 16 +- 2d_fourie/test_voronoi.py | 46 +++++ 7 files changed, 459 insertions(+), 337 deletions(-) create mode 100644 2d_fourie/analysis.py create mode 100644 2d_fourie/extractors.py create mode 100644 2d_fourie/plot.py create mode 100644 2d_fourie/test_voronoi.py diff --git a/2d_fourie/analysis.py b/2d_fourie/analysis.py new file mode 100644 index 0000000..208a18a --- /dev/null +++ b/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() diff --git a/2d_fourie/extractors.py b/2d_fourie/extractors.py new file mode 100644 index 0000000..f296d71 --- /dev/null +++ b/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() diff --git a/2d_fourie/lattices.py b/2d_fourie/lattices.py index 8074e0b..b3b17cf 100644 --- a/2d_fourie/lattices.py +++ b/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 @@ -28,10 +29,10 @@ 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,Y = np.meshgrid(x, y) - return [(X,Y)] + x = np.arange(-3, 3) * 0.2 + y = np.arange(-3, 3) * 0.2 + X, Y = np.meshgrid(x, y) + return [(X, Y)] class VO2_Lattice(Lattice): @@ -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,19 +104,28 @@ class VO2_Lattice(Lattice): return inplace_pos_x, inplace_pos_y def reci_rutile(self): - 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() - + x = np.arange(-20, 21) + 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 + x, y = self.reci_rutile() + return x + 0.1083, y + 0.1719 def reci_mono_2(self): - x, y = self.reci_rutile() - return x - 0.1083, y + 0.1719 + x, y = self.reci_rutile() + return x - 0.1083, y + 0.1719 def reci(self): - return [self.reci_rutile(), self.reci_mono(), self.reci_mono_2()] \ No newline at end of file + 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] diff --git a/2d_fourie/main.py b/2d_fourie/main.py index 54ede8a..2c18e14 100644 --- a/2d_fourie/main.py +++ b/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" ) @@ -101,11 +63,11 @@ def rotate(x, y, angle): def test_square(): LEN = 40 - #lat = SCC_Lattice(LEN, LEN) + # lat = SCC_Lattice(LEN, LEN) lat = VO2_Lattice(LEN, LEN) plot = Plotter(lat) pos_x, pos_y = lat.get_from_mask(np.zeros((40, 40))) - #pos_x, pos_y = rotate(pos_x, pos_y, 30) + # pos_x, pos_y = rotate(pos_x, pos_y, 30) si = SpinImage(pos_x, pos_y) fig, axs = plt.subplots(2, 2) si.pad_it_square(10) @@ -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: - plt.scatter(percentage, o/o[0]) - plt.plot([0,1], [o[0], o[-1]]) - plt.show() + 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) if __name__ == "__main__": # test_square() - # test_mixed() - random() -# def test_lattice(): -# lat = VO2_Lattice(10, 10) -# maske = np.zeros((10, 10), dtype=bool) -# x, y = lat.get_from_mask(maske) -# -# plt.scatter(x, y) -# maske = np.invert(maske) -# 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() + test_mixed() + plt.show() + # random() + np.random.seed(1234) + for i in np.random.randint(0, 10000, 2): + random(i) + ising(i) + + plt.show() diff --git a/2d_fourie/plot.py b/2d_fourie/plot.py new file mode 100644 index 0000000..5b685d0 --- /dev/null +++ b/2d_fourie/plot.py @@ -0,0 +1,5 @@ +import matplotlib.pyplot as plt +import numpy as np + +if __name__ == "__main__": + pass diff --git a/2d_fourie/spin_image.py b/2d_fourie/spin_image.py index 39c678a..6980f08 100644 --- a/2d_fourie/spin_image.py +++ b/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) diff --git a/2d_fourie/test_voronoi.py b/2d_fourie/test_voronoi.py new file mode 100644 index 0000000..9ba2806 --- /dev/null +++ b/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()