diff --git a/2d_fourie/cache.py b/2d_fourie/cache.py index 389b6d4..be741e3 100644 --- a/2d_fourie/cache.py +++ b/2d_fourie/cache.py @@ -7,12 +7,13 @@ def timeit(func): @wraps(func) def timeit_wrapper(*args, **kwargs): start_time = time.perf_counter() + print(f"Start Function {func.__name__}:") result = func(*args, **kwargs) end_time = time.perf_counter() total_time = end_time - start_time # first item in the args, ie `args[0]` is `self` print( - f'Function {func.__name__}{args} {kwargs} Took {total_time:.4f} seconds') + f'Function {func.__name__} Took {total_time:.4f} seconds') return result return timeit_wrapper diff --git a/2d_fourie/lattices.py b/2d_fourie/lattices.py index 2d09b71..d25e21b 100644 --- a/2d_fourie/lattices.py +++ b/2d_fourie/lattices.py @@ -16,6 +16,9 @@ class Lattice: def get_from_mask(self, maske): pass + def get_both(self): + pass + def reci(self): pass @@ -30,6 +33,9 @@ class SCC_Lattice(Lattice): def get_from_mask(self, maske): return self.X, self.Y + def get_both(self): + return [(self.X, self.Y), (self.X, self.Y)] + def reci(self): x = np.arange(-3, 4) * 0.2 y = np.arange(-3, 4) * 0.2 @@ -71,8 +77,6 @@ class VO2_Lattice(Lattice): offset_a_r = res * int(offset_a_r/res) offset_c_r = res * int(offset_c_r/res) - offset_a_r = 0.5 - offset_c_r = 0.5 print(offset_a_r, offset_c_r) self.atom_x_mono = offset_a_r + X * \ @@ -116,9 +120,11 @@ class VO2_Lattice(Lattice): inplace_pos_y[mask] = self.atom_y_mono[mask] return inplace_pos_x, inplace_pos_y + def get_both(self): + return [(self.atom_x_rut, self.atom_y_rut), (self.atom_x_mono, self.atom_y_mono)] + def reci_rutile(self): num = 20 - #num = 2 x = np.arange(-num, num + 1) y = np.arange(-num, num + 1) X, Y = np.meshgrid(x, y) diff --git a/2d_fourie/main.py b/2d_fourie/main.py index 5e46d4c..85caba1 100644 --- a/2d_fourie/main.py +++ b/2d_fourie/main.py @@ -1,5 +1,5 @@ from lattices import SCC_Lattice, VO2_Lattice -from spin_image import SpinImage +from spin_image import SpinImage, SpinImage_Two_Phase import numpy as np import matplotlib.pyplot as plt import matplotlib.patches as patches @@ -7,6 +7,7 @@ import matplotlib import tqdm from extractors import Rect_Evaluator, Voronoi_Evaluator, Image_Wrapper from cache import timeit +from scipy import signal class Plotter: @@ -38,7 +39,7 @@ class Plotter: intens, extent=(np.min(freqx), np.max(freqx), np.min(freqy), np.max(freqy)), - norm=matplotlib.colors.LogNorm(vmin=10), + norm=matplotlib.colors.LogNorm(vmin=1e-12), cmap="viridis", origin="lower" ) @@ -55,6 +56,14 @@ class Plotter: plt.colorbar(t, ax=ax_lin) +def plot_periodic(array): + plt.figure() + kernel = np.ones((10, 10)) + #array = signal.convolve2d(array, kernel, boundary='symm', mode='same') + tiled = np.tile(array, (2, 2)) + plt.imshow(tiled) + + def rotate(x, y, angle): radian = angle / 180 * 2 * np.pi return np.cos(radian) * x - np.sin(radian) * y,\ @@ -67,14 +76,18 @@ def test_square(): # lat = VO2_Lattice(LEN, LEN) plot = Plotter(lat) pos_x, pos_y = lat.get_from_mask(np.zeros((40, 40))) + [(x, y), (x1, y1)] = lat.get_both() # pos_x, pos_y = rotate(pos_x, pos_y, 30) - si = SpinImage(pos_x, pos_y) + # si = SpinImage(pos_x, pos_y) + si = SpinImage_Two_Phase(x, y, x1, y1) + si.apply_mask(np.zeros((20, 20))) fig, axs = plt.subplots(2, 2) - si.plot(axs[0, 0], 1) + si.plot(axs[0, 0]) # si.gaussian(LEN) # si.blur(3) - si.pad_it_square(size=int((5*LEN)/si.resolution)) - si.plot(axs[0, 1], 1) + # si.pad_it_square(additional_pad=1000) + si.gaussian(20) + si.plot(axs[0, 1]) plt.pause(0.1) @@ -139,10 +152,14 @@ def test_mixed(): int((LEN * 4.554) / si_mixed.resolution)+shift_y) fx, fy, intens_mixed = si_mixed.fft() + plot_periodic(si_rutile.img) + plot_periodic(si_mono.img) + plot_periodic(si_mixed.img) + fig, axs = plt.subplots(3, 3) - si_rutile.plot(axs[0, 0], 1) - si_mono.plot(axs[0, 2], 1) - si_mixed.plot(axs[0, 1], 1) + si_rutile.plot(axs[0, 0]) + si_mono.plot(axs[0, 2]) + si_mixed.plot(axs[0, 1]) plot.plot(freqx=fx, freqy=fy, intens=intens_rutile, ax_log=axs[1, 0], ax_lin=axs[2, 0], vmax=10e7) plot.plot(freqx=fx, freqy=fy, intens=intens_mono, @@ -295,9 +312,37 @@ def ising(seed): out_2=out_voro[1], out_3=out_voro[2], out_4=out_voro[3]) +def test_me(): + LEN = 20 + lat = VO2_Lattice(LEN, LEN) + maske = np.zeros((LEN, LEN), dtype=bool) + pos_x, pos_y = lat.get_from_mask(maske) + si = SpinImage(pos_x, pos_y) + + fig, axs = plt.subplots(1, 3) + kernel = np.ones((20, 20)) + array = signal.convolve2d(si.img, kernel, boundary='symm', mode='same') + axs[0].imshow(array) + + [(x, y), (x1, y1)] = lat.get_both() + si = SpinImage_Two_Phase(x, y, x1, y1) + si.apply_mask(maske) + + array = signal.convolve2d(si.img, kernel, boundary='symm', mode='same') + # axs[1].imshow(array) + axs[1].imshow(si.img) + + np.invert(maske) + si.apply_mask(maske) + array = signal.convolve2d(si.img, kernel, boundary='symm', mode='same') + # axs[1].imshow(array) + axs[2].imshow(si.img) + + if __name__ == "__main__": - # test_square() - test_mixed() + # test_me() + test_square() + # test_mixed() plt.show() # random() # np.random.seed(1234) diff --git a/2d_fourie/spin_image.py b/2d_fourie/spin_image.py index d81ca58..bb72597 100644 --- a/2d_fourie/spin_image.py +++ b/2d_fourie/spin_image.py @@ -2,6 +2,183 @@ import numpy as np import scipy import scipy.fftpack as sfft import scipy.signal +import tqdm +from cache import timeit + + +class SpinImage_Two_Phase: + resolution = 0.05 + offset = 40 + + @timeit + def __init__(self, x_pos_low, y_pos_low, x_pos_high, y_pos_high): + assert x_pos_low.shape == y_pos_low.shape + assert x_pos_high.shape == y_pos_high.shape + assert x_pos_low.shape == x_pos_high.shape + + offset_shift = self.offset * self.resolution + + x_pos_low = x_pos_low - np.min(x_pos_low) + offset_shift + y_pos_low = y_pos_low - np.min(y_pos_low) + offset_shift + x_pos_high = x_pos_high - np.min(x_pos_high) + offset_shift + y_pos_high = y_pos_high - np.min(y_pos_high) + offset_shift + max_len_x = np.maximum(np.max(x_pos_low), np.max(x_pos_high)) + max_len_y = np.maximum(np.max(y_pos_low), np.max(y_pos_high)) + self.length_x = max_len_x + (self.offset + 1) * self.resolution + self.length_y = max_len_y + (self.offset + 1) * self.resolution + self.x_ind = np.arange(0, self.length_x, self.resolution) # angstrom + self.y_ind = np.arange(0, self.length_y, self.resolution) # angstrom + self.x_low, self.y_low = self._stuff(x_pos_low, y_pos_low) + self.x_high, self.y_high = self._stuff(x_pos_high, y_pos_high) + + X, Y = np.meshgrid(self.x_ind, self.y_ind, indexing="ij") + print(X.shape, self.x_low.flatten().max(), self.y_low.flatten().max()) + sigma = .1 + self.low_list = [] + for x, y, x_ind, y_ind in zip( + x_pos_low.flatten(), + y_pos_low.flatten(), + self.x_low.flatten(), + self.y_low.flatten(), + ): + xl, yl, xu, yu = self.clean_bounds( + x_ind - self.offset, + y_ind - self.offset, + x_ind + self.offset, + y_ind + self.offset, + X.shape, + ) + self.low_list.append(np.exp(-0.5 * ((X[xl:xu, yl:yu] - x) ** 2 + + (Y[xl:xu, yl:yu] - y) ** 2) / sigma**2)) + self.low_list = np.array( + self.low_list, dtype=object) + #print("BEFOR:", self.low_list.shape) + #self.low_list = self.low_list.reshape((*self.x_high.shape,80,80)) + #print("HIER:", self.low_list.shape) + self.high_list = [] + for x, y, x_ind, y_ind in zip( + x_pos_high.flatten(), + y_pos_high.flatten(), + self.x_high.flatten(), + self.y_high.flatten(), + ): + xl, yl, xu, yu = self.clean_bounds( + x_ind - self.offset, + y_ind - self.offset, + x_ind + self.offset, + y_ind + self.offset, + X.shape, + ) + self.high_list.append(np.exp(-0.5 * ((X[xl:xu, yl:yu] - x) ** 2 + + (Y[xl:xu, yl:yu] - y) ** 2) / sigma**2)) + self.high_list = np.array(self.high_list, dtype=object) + + def clean_bounds(self, xl, yl, xu, yu, shape): + if xl < 0: + xl = 0 + if yl < 0: + yl = 0 + + if xu > shape[0]: + xu = shape[0] + if yu > shape[1]: + yu = shape[1] + + return xl, yl, xu, yu + + def _stuff(self, pos_x, pos_y): + xind = np.searchsorted(self.x_ind, pos_x).astype(int) + yind = np.searchsorted(self.y_ind, pos_y).astype(int) + return xind, yind + + @timeit + def apply_mask(self, maske): + mask = np.empty_like(self.x_high, dtype=bool) + mask[0::2, 0::2] = maske + mask[1::2, 0::2] = maske + mask[0::2, 1::2] = maske + mask[1::2, 1::2] = maske + mask = mask.flatten() + self.img = np.zeros((self.x_ind.size, self.y_ind.size)) + print(self.img.shape) + for x, y, dat in zip( + self.x_high.flatten()[mask], + self.y_high.flatten()[mask], + self.high_list[mask], + ): + xl, yl, xu, yu = self.clean_bounds( + x - self.offset, + y - self.offset, + x + self.offset, + y + self.offset, + self.img.shape, + ) + self.img[x - self.offset: x + self.offset, + y - self.offset: y + self.offset] = dat + self.img[x, y] += 1 + mask = np.invert(mask) + for x, y, dat in zip( + self.x_high.flatten()[mask], + self.y_high.flatten()[mask], + self.high_list[mask], + ): + xl, yl, xu, yu = self.clean_bounds( + x - self.offset, + y - self.offset, + x + self.offset, + y + self.offset, + self.img.shape, + ) + if self.img[xl:xu, yl:yu].shape == dat.shape: + self.img[xl:xu, yl:yu] = dat + else: + self.img[x, y] += 1 + + @timeit + def fft(self): + Z_fft = sfft.fft2(self.img) + Z_shift = sfft.fftshift(Z_fft) + fft_freqx = sfft.fftfreq(self.img.shape[0], self.resolution) + fft_freqy = sfft.fftfreq(self.img.shape[1], self.resolution) + fft_freqx_clean = sfft.fftshift(fft_freqx) + fft_freqy_clean = sfft.fftshift(fft_freqy) + return fft_freqx_clean, fft_freqy_clean, np.abs(Z_shift) ** 2 + + @timeit + def plot(self, ax, scale=None): + if scale is None: + ax.imshow(self.img) + else: + quad = np.ones((int(scale / self.resolution), + int(scale / self.resolution))) + img = scipy.signal.convolve2d(self.img, quad) + ax.imshow(img) + + @timeit + 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 + 1 + if size is not None: + xx = np.maximum(xx, size) + yy = xx + + a = (xx - h) // 2 + aa = xx - a - h + + b = (yy - w) // 2 + bb = yy - b - w + + self.img = np.pad(self.img, pad_width=( + (a, aa), (b, bb)), mode="constant") + + def gaussian(self, sigma): + x = np.arange(-self.length_x / 2, self.length_x / 2, self.resolution) + y = np.arange(-self.length_y / 2, self.length_y / 2, self.resolution) + 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))) + self.img = np.multiply(self.img, z.T) class SpinImage: @@ -23,6 +200,15 @@ class SpinImage: img[xind, yind] = 1 return img + def image_from_pos_with_gauss(self, pos_x, pos_y, sigma=1): + x_ind = np.arange(0, self.length_x, self.resolution) # angstrom + y_ind = np.arange(0, self.length_y, self.resolution) # angstrom + img = np.zeros((x_ind.size, y_ind.size)) + X, Y = np.meshgrid(y_ind, x_ind) + for px, py in tqdm.tqdm(zip(pos_x.flatten(), pos_y.flatten())): + img += np.exp(-0.5 * ((X - px) ** 2 + (Y - py) ** 2) / sigma**2) + return img + def fft(self): Z_fft = sfft.fft2(self.img) Z_shift = sfft.fftshift(Z_fft) @@ -69,27 +255,19 @@ class SpinImage: (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 = 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))) - ) + 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) - y = np.arange(-self.length_y / 2, - self.length_y / 2, self.resolution) + x = np.arange(-self.length_x / 2, self.length_x / 2, self.resolution) + y = np.arange(-self.length_y / 2, self.length_y / 2, self.resolution) 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))) - ) + z = 1 / (2 * np.pi * sigma * sigma) * \ + np.exp(-(X**2 / (2 * sigma**2) + Y**2 / (2 * sigma**2))) self.img = np.multiply(self.img, z.T) def plot(self, ax, scale=None): @@ -97,7 +275,7 @@ class SpinImage: ax.imshow(self.img) else: quad = np.ones((int(scale / self.resolution), - int(scale / self.resolution))) + int(scale / self.resolution))) img = scipy.signal.convolve2d(self.img, quad) ax.imshow(img)