diff --git a/2d_fourie/main.py b/2d_fourie/main.py index 85caba1..9da3b09 100644 --- a/2d_fourie/main.py +++ b/2d_fourie/main.py @@ -121,45 +121,32 @@ def helper(intens, fx, fy, voro, rect): def test_mixed(): - shift_x = -(5 + 3*28) - shift_y = -43 - LEN = 40 + fig, axs = plt.subplots(3, 3) + LEN = 80 lat = VO2_Lattice(LEN, LEN) plot = Plotter(lat) - - pos_x, pos_y = lat.get_from_mask(np.zeros((LEN, LEN))) - si_mono = SpinImage(pos_x, pos_y) - si_mono.pad_it(int((LEN * 5.712) / si_mono.resolution)+shift_x, - int((LEN * 4.554) / si_mono.resolution)+shift_y) - fx, fy, intens_mono = si_mono.fft() - - pos_x, pos_y = lat.get_from_mask(np.ones((LEN, LEN))) - si_rutile = SpinImage(pos_x, pos_y) - print(int((LEN * 5.712) / si_rutile.resolution)+shift_x, - int((LEN * 4.554) / si_rutile.resolution)+shift_y, - si_rutile.img.shape) - si_rutile.pad_it(int((LEN * 5.712) / si_rutile.resolution)+shift_x, - int((LEN * 4.554) / si_rutile.resolution)+shift_y) - fx, fy, intens_rutile = si_rutile.fft() - - mask_misk = np.ones((LEN, LEN)) + [(x, y), (x1, y1)] = lat.get_both() + si = SpinImage_Two_Phase(x, y, x1, y1) + mask_misk = np.ones((LEN, LEN), dtype=bool) ind = np.arange(mask_misk.size) np.random.shuffle(ind) mask_misk[np.unravel_index(ind[:800], (LEN, LEN))] = False - pos_x, pos_y = lat.get_from_mask(mask_misk) - si_mixed = SpinImage(pos_x, pos_y) - si_mixed.pad_it(int((LEN * 5.712) / si_mixed.resolution)+shift_x, - 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) + si.apply_mask(np.zeros((LEN, LEN), dtype=bool)) + si.gaussian(20) + fx, fy, intens_mono = si.fft() + si.plot(axs[0, 2]) + + si.apply_mask(np.ones((LEN, LEN))) + si.gaussian(20) + fx, fy, intens_rutile = si.fft() + si.plot(axs[0, 0]) + + si.apply_mask(mask_misk) + si.gaussian(20) + fx, fy, intens_mixed = si.fft() + si.plot(axs[0, 1]) - fig, axs = plt.subplots(3, 3) - 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, @@ -167,30 +154,13 @@ def test_mixed(): plot.plot(freqx=fx, freqy=fy, intens=intens_mixed, ax_log=axs[1, 1], ax_lin=axs[2, 1], vmax=10e7) - index = np.abs(fx).argmin() - print(index, "/", fx.shape) - fig, axs = plt.subplots(1, 3) - axs[0].plot(intens_rutile[index, :]) - axs[0].plot(intens_rutile[index-1, :]) - axs[0].plot(intens_rutile[index+1, :]) - axs[1].plot(intens_mixed[index, :]) - axs[2].plot(intens_mono[index, :]) - for a in axs: - a.set_yscale("log") + # Plotting cuts - index = np.abs(fy).argmin() - print(index, "/", fy.shape) - fig, axs = plt.subplots(1, 3) - axs[0].plot(intens_rutile[:, index]) - axs[1].plot(intens_mixed[:, index]) - axs[2].plot(intens_mono[:, index]) - for a in axs: - a.set_yscale("log") def random(seed): np.random.seed(seed) - LEN = 40 + LEN = 80 lat = VO2_Lattice(LEN, LEN) maske = np.zeros((LEN, LEN)) ind = np.arange(LEN * LEN) @@ -319,7 +289,7 @@ def test_me(): pos_x, pos_y = lat.get_from_mask(maske) si = SpinImage(pos_x, pos_y) - fig, axs = plt.subplots(1, 3) + fig, axs = plt.subplots(1, 4) kernel = np.ones((20, 20)) array = signal.convolve2d(si.img, kernel, boundary='symm', mode='same') axs[0].imshow(array) @@ -328,21 +298,19 @@ def test_me(): 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) + tmp = si.img + maske = 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) + axs[3].imshow(si.img+tmp) + if __name__ == "__main__": # test_me() - test_square() - # test_mixed() + # 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 bb72597..d2d26b4 100644 --- a/2d_fourie/spin_image.py +++ b/2d_fourie/spin_image.py @@ -10,6 +10,31 @@ class SpinImage_Two_Phase: resolution = 0.05 offset = 40 + def make_list(self, x_pos, y_pos, x_inds, y_inds): + sigma = .1 + x_ind = np.arange(0, self.length_x, self.resolution) + y_ind = np.arange(0, self.length_y, self.resolution) + X, Y = np.meshgrid(x_ind, y_ind, indexing="ij") + out_list = [] + for x, y, x_ind, y_ind in zip( + x_pos.flatten(), + y_pos.flatten(), + x_inds.flatten(), + y_inds.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, + ) + out_list.append(np.exp(-0.5 * ((X[xl:xu, yl:yu] - x) ** 2 + + (Y[xl:xu, yl:yu] - y) ** 2) / sigma**2)) + #out_list.append(np.ones_like(X[xl:xu, yl:yu])) + out_list = np.array(out_list) + return out_list + @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 @@ -18,60 +43,28 @@ class SpinImage_Two_Phase: 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 + zero_shift_x = np.minimum(np.min(x_pos_low), np.min(x_pos_high)) + zero_shift_y = np.minimum(np.min(y_pos_low), np.min(y_pos_high)) + + x_pos_low = x_pos_low - zero_shift_x + offset_shift + y_pos_low = y_pos_low - zero_shift_y + offset_shift + + x_pos_high = x_pos_high - zero_shift_x + offset_shift + y_pos_high = y_pos_high - zero_shift_y + 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) + self.x_index_low, self.y_index_low = self._stuff(x_pos_low, y_pos_low) + self.x_index_high, self.y_index_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) + self.buffer_low = self.make_list( + x_pos_low, y_pos_low, self.x_index_low, self.y_index_low) + self.buffer_high = self.make_list( + x_pos_high, y_pos_high, self.x_index_high, self.y_index_high) def clean_bounds(self, xl, yl, xu, yu, shape): if xl < 0: @@ -87,24 +80,17 @@ class SpinImage_Two_Phase: 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) + x_ind = np.arange(0, self.length_x, self.resolution) + y_ind = np.arange(0, self.length_y, self.resolution) + xind = np.searchsorted(x_ind, pos_x).astype(int) + yind = np.searchsorted(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) + def _apply_mask(self, x, y, buffer, mask): for x, y, dat in zip( - self.x_high.flatten()[mask], - self.y_high.flatten()[mask], - self.high_list[mask], + x.flatten()[mask], + y.flatten()[mask], + buffer[mask], ): xl, yl, xu, yu = self.clean_bounds( x - self.offset, @@ -114,25 +100,24 @@ class SpinImage_Two_Phase: self.img.shape, ) self.img[x - self.offset: x + self.offset, - y - self.offset: y + self.offset] = dat - self.img[x, y] += 1 + y - self.offset: y + self.offset] += dat + + @timeit + def apply_mask(self, maske): + mask = np.empty_like(self.x_index_high, dtype=bool) + print(maske) + 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( + (int(self.length_x/self.resolution), int(self.length_y/self.resolution))) + self._apply_mask(self.x_index_low, self.y_index_low, + self.buffer_low, mask) 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 + self._apply_mask(self.x_index_high, self.y_index_high, + self.buffer_high, mask) @timeit def fft(self): @@ -173,9 +158,11 @@ class SpinImage_Two_Phase: (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 = np.linspace(0, self.length_x, self.img.shape[0]) + y = np.linspace(0, self.length_y, self.img.shape[1]) X, Y = np.meshgrid(x, y) + X = X - (self.length_x / 2.) + Y = Y - (self.length_y / 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)