Jump Probability with Greyscale Images

checkers_and_stripes.py

@@ -0,0 +1,42 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue May 12 16:09:18 2020
+
+@author: Tom
+"""
+
+import pytrax as pt
+import numpy as np
+import matplotlib.pyplot as plt
+plt.close('all')
+if __name__ == '__main__':
+ n = 5
+ im_chess = np.ones([10*n, 10*n])
+ g = 0.25
+ for i in range(10*n):
+ for j in range(10*n):
+ if (np.floor(i/n) % 2) != (np.floor(j/n) % 2):
+ im_chess[i, j] = g
+ plt.figure()
+ plt.imshow(im_chess)
+ rw_c = pt.RandomWalk(im_chess)
+ rw_c.run(nt=10000, nw=10000, num_proc=1, same_start=True)
+ rw_c.calc_msd()
+ rw_c.plot_msd()
+ rw_c.plot_walk_2d()
+ rw_c.plot_walk_2d(w_id=np.argmin(rw_c.sq_disp[-1, :]), data='t')
+ rw_c.plot_walk_2d(w_id=np.argmax(rw_c.sq_disp[-1, :]), data='t')
+
+ im = np.ones([10*n, 10*n])
+ for i in range(10*n):
+ if (np.floor(i/n) % 2 == 0):
+ im[i, :] = g
+ plt.figure()
+ plt.imshow(im)
+ rw = pt.RandomWalk(im)
+ rw.run(nt=10000, nw=10000, stride=10, num_proc=1, same_start=True)
+ rw.calc_msd()
+ rw.plot_msd()
+ rw.plot_walk_2d()
+ rw.plot_walk_2d(w_id=np.argmin(rw.sq_disp[-1, :]), data='t')
+ rw.plot_walk_2d(w_id=np.argmax(rw.sq_disp[-1, :]), data='t')

dots.py

@@ -0,0 +1,114 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jun 3 11:07:23 2020
+
+@author: tom
+"""
+
+if __name__ == '__main__':
+ import numpy as np
+ import matplotlib.pyplot as plt
+ import scipy.ndimage as spim
+ import pytrax as pt
+ from scipy.interpolate import NearestNDInterpolator
+
+ dots = np.zeros([1000, 1000])
+ for i in range(10):
+ for j in range(10):
+ adjx = np.random.choice(np.arange(-25, 25, 1, int), 1)[0]
+ adjy = np.random.choice(np.arange(-25, 25, 1, int), 1)[0]
+ dots[(100*i)+50+adjx, (100*j)+50+adjy] = 1
+
+ dt = spim.distance_transform_edt(dots)
+ strel = spim.generate_binary_structure(2, 1)
+ big_dots = spim.morphology.binary_dilation(dots, strel, 10)
+ plt.figure()
+ plt.imshow(big_dots)
+ dt = spim.distance_transform_edt(1-big_dots)
+ dt = dt/dt.max()
+ dt = 1 - dt
+ plt.figure()
+ plt.imshow(dt)
+ plt.figure()
+ plt.hist(dt.flatten())
+
+ def grey_scaler(im, lower, upper, exponent):
+ im_scaled = im.copy().astype(float)
+ # below lower is pore
+ im_scaled[im < lower] = lower
+ # above upper is nmc
+ im_scaled[im > upper] = upper
+ im_scaled -= lower
+ # normalize scale
+ im_scaled /= (upper-lower)
+ # invert image
+ im_scaled = (1-im_scaled)
+ # apply exponent to grey values - higher exponent means slower transport in regions near nmc intensity
+ im_scaled = im_scaled**exponent
+ return im_scaled
+
+ def process(im, thresh):
+ tmp = im < thresh
+ tmp[:10, :] = True
+ tmp[-10:, :] = True
+ lab, N = spim.label(tmp)
+ return lab == 1
+
+ def process_scaled(im, thresh, trange):
+ tmp = grey_scaler(im, thresh, thresh+trange, 1.0)
+ tmp[:10, :] = 1.0
+ tmp[-10:, :] = 1.0
+ tmp_bin = tmp > 0.0
+ lab, N = spim.label(tmp_bin)
+ tmp[lab > 1] = 0.0
+ return tmp
+
+ def tortuosity(image):
+ rw = pt.RandomWalk(image)
+ rw.run(nt=100000, nw=10000, stride=10, num_proc=10)
+ rw.calc_msd()
+ rw.plot_msd()
+ return rw
+
+ def interpolated_sq_disp(rw):
+ rw.colour_sq_disp()
+ im = rw.im_sq_disp
+ x_len, y_len = im.shape
+ points = np.argwhere(im > 1)
+ colours = im[points[:, 0], points[:, 1]]
+ myInterpolator = NearestNDInterpolator(points, colours)
+ grid_x, grid_y = np.mgrid[0:x_len:np.complex(x_len, 0),
+ 0:y_len:np.complex(y_len, 0)]
+ arr = np.log(myInterpolator(grid_x, grid_y).astype(float))
+ arr[im == 0] = np.nan
+ plt.figure()
+ plt.imshow(arr)
+ plt.colorbar()
+
+ plt.figure()
+ plt.imshow(dt)
+ # thresh = np.linspace(0.525, 0.55, 3)
+ thresh = [0.6]
+ rws = []
+ ims = []
+ for i in range(len(thresh)):
+ plt.figure()
+ tmp = process(dt, thresh[i])
+ plt.imshow(tmp)
+ plt.title('Threshold: '+str(np.around(thresh[i], 3))+' Porosity: '+str(np.around(np.sum(tmp)/np.size(tmp), 3)))
+ rws.append(tortuosity(tmp))
+ ims.append(tmp)
+ tau_0 = [r.data['axis_0_tau'] for r in rws]
+ # grey_rws = []
+ # grey_ims = []
+ # for i in range(len(thresh)):
+ # plt.figure()
+ # tmp = process_scaled(dt, thresh[i], 0.05)
+ # plt.imshow(tmp)
+ # plt.title('Threshold: '+str(np.around(thresh[i], 3))+' Porosity: '+str(np.around(np.sum(tmp)/np.size(tmp), 3)))
+ # grey_rws.append(tortuosity(tmp))
+ # grey_ims.append(tmp)
+ # grey_tau_0 = [r.data['axis_0_tau'] for r in grey_rws]
+ plt.figure()
+ plt.plot(tau_0)
+ # plt.plot(grey_tau_0)

frequency.py

@@ -0,0 +1,85 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Oct 10 08:34:14 2019
+
+@author: Tom
+"""
+
+import porespy as ps
+import pytrax as pt
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.ndimage as spim
+import time
+
+plt.close('all')
+def main():
+ im = ps.generators.blobs(shape=[1000, 1000], blobiness=3, porosity=0.5)
+ im = ps.filters.fill_blind_pores(im).astype(int)
+ dt = spim.distance_transform_edt(im)
+ grey = dt.copy()/dt.max()
+ grey = np.pad(grey, 1, mode='constant', constant_values=0)
+ # Number of time steps and walkers
+ num_t = 100000
+ num_w = 1000
+ stride = 1
+
+ rw = pt.RandomWalk(grey, seed=False)
+ rw.run(num_t, num_w, same_start=False, stride=stride, num_proc=12)
+ # Plot mean square displacement
+ rw.plot_msd()
+ rw.plot_walk_2d()
+
+ print('Calculating hit frequency')
+ coords = rw.real_coords
+ freq = np.zeros_like(grey)
+ if len(im.shape) == 2:
+ x_last = coords[0, :, 0].fill(-1)
+ y_last = coords[0, :, 1].fill(-1)
+ for t in range(coords.shape[0]):
+ x = coords[t, :, 0]
+ y = coords[t, :, 1]
+ same_x = x == x_last
+ same_y = y == y_last
+ same_xy = same_x * same_y
+ freq[x[~same_xy], y[~same_xy]] += 1
+ x_last = x
+ y_last = y
+ else:
+ x_last = coords[0, :, 0].fill(-1)
+ y_last = coords[0, :, 1].fill(-1)
+ z_last = coords[0, :, 2].fill(-1)
+ for t in range(coords.shape[0]):
+ x = coords[t, :, 0]
+ y = coords[t, :, 1]
+ z = coords[t, :, 2]
+ same_x = x == x_last
+ same_y = y == y_last
+ same_z = z == z_last
+ same_xyz = same_x * same_y * same_z
+ freq[x[~same_xyz], y[~same_xyz], z[~same_xyz]] += 1
+ x_last = x
+ y_last = y
+ z_last = z
+ return grey, freq
+
+if __name__ == '__main__':
+ st = time.time()
+ grey, freq = main()
+ some_hits = freq[freq > 0]
+ frange = np.unique(some_hits)
+ plt.figure()
+ plt.hist(some_hits, bins=int(frange.max()-frange.min()))
+ log_freq = np.log(freq)
+ log_freq[freq == 0] = np.nan
+ freq[freq == 0] = np.nan
+ print(frange.min(), frange.max())
+ plt.figure()
+ plt.imshow(grey > 0, cmap='gist_gray')
+ plt.imshow(freq)
+ plt.colorbar()
+ plt.figure()
+ plt.imshow(grey > 0, cmap='gist_gray')
+ plt.imshow(log_freq)
+ plt.colorbar()
+ print('Sim time', time.time() - st)

grey.py

@@ -0,0 +1,32 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Oct 8 17:07:13 2019
+
+@author: Tom
+"""
+
+import porespy as ps
+import pytrax as pt
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.ndimage as spim
+if __name__ == '__main__':
+ im = ps.generators.blobs(shape=[200, 200], porosity=0.5).astype(int)
+ plt.figure()
+ plt.imshow(im)
+ dt = spim.distance_transform_edt(im)
+ grey = dt.copy()/dt.max()
+ # Number of time steps and walkers
+ num_t = 100000
+ num_w = None
+ stride = 10
+
+ for case in [im]:
+ rw = pt.RandomWalk(case, seed=False)
+ rw.run(num_t, num_w, same_start=False, stride=stride, num_proc=10)
+ # Plot mean square displacement
+ rw.plot_msd()
+ # rw.plot_walk_2d()
+ rw.colour_sq_disp()
+ plt.figure()
+ plt.imshow(rw.im_sq_disp)

multi_grey.py

@@ -0,0 +1,30 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 24 10:51:10 2020
+
+@author: Tom
+"""
+
+import pytrax as pt
+import numpy as np
+import matplotlib.pyplot as plt
+plt.close('all')
+
+if __name__ == '__main__':
+ im = np.ones([1000, 1000])
+ # Number of time steps and walkers
+ num_t = 10000
+ num_w = 800
+ stride = 1
+ grey_vals = np.logspace(-1, 1, 10)
+ tau = []
+ for grey_val in grey_vals:
+ grey = im.copy()
+ grey = grey.astype(float)
+ grey[grey == 1.0] = grey_val
+ rw = pt.RandomWalk(grey, seed=False)
+ rw.run(num_t, num_w, same_start=False, stride=stride, num_proc=1)
+ rw.plot_msd()
+ tau.append(rw.data['axis_0_tau'])
+ plt.figure()
+ plt.loglog(grey_vals, tau)