small changes

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, np.int), 1)[0]
+ adjy = np.random.choice(np.arange(-25, 25, 1, np.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.525]
+ 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)

grey.py

@@ -11,17 +11,22 @@
 import matplotlib.pyplot as plt
 import scipy.ndimage as spim
 if __name__ == '__main__':
- im = ps.generators.blobs(shape=[1000, 1000], porosity=0.7).astype(np.int)
+ im = ps.generators.blobs(shape=[200, 200], porosity=0.5).astype(np.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 = 10000
- num_w = 800
- stride = 1
+ num_t = 100000
+ num_w = None
+ stride = 10
 
- for case in [im, grey]:
+ 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.plot_walk_2d()
+ rw.colour_sq_disp()
+ plt.figure()
+ plt.imshow(rw.im_sq_disp)

pytrax/__RandomWalk__.py

@@ -192,7 +192,7 @@ def check_edge(self, walkers, axis, move, real):
 
  return move, move_real, real
 
- def _get_starts(self, same_start=False):
+ def _get_starts(self, same_start=False, max_iter=100):
  r'''
  Start walkers in the pore space at random location
  same_start starts all the walkers at the same spot if True and at
@@ -203,7 +203,7 @@ def _get_starts(self, same_start=False):
  determines whether to start all the walkers at the same coordinate
  '''
  if not same_start:
- walkers = self._rand_start(self.im, num=self.nw)
+ walkers = self._rand_start(self.im, num=self.nw) 
  else:
  w = self._rand_start(self.im, num=1).flatten()
  walkers = np.tile(w, (self.nw, 1))
@@ -255,6 +255,7 @@ def _run_walk(self, walkers):
  real = np.ones_like(walkers)
 # real_coords = np.ndarray([self.nt, nw, self.dim], dtype=int)
  real_coords = []
+ real_coords.append(wr.copy())
  for t in range(self.nt):
  # Random velocity update
  # Randomly select an axis to move along for each walker
@@ -301,7 +302,8 @@ def run(self, nt=1000, nw=1, same_start=False, stride=1, num_proc=1):
  nt: int (default = 1000)
  the number of timesteps to run the simulation for
  nw: int (default = 1)
- he vector of the next move to be made by the walker
+ the number of walkers to run in the simulation - None will start 1
+ walker in every non-zero pixel of the image
  same_start: bool
  determines whether to start all the walkers at the same coordinate
  stride: int
@@ -312,11 +314,16 @@ def run(self, nt=1000, nw=1, same_start=False, stride=1, num_proc=1):
  multiprocessing.
  '''
  self.nt = int(nt)
- self.nw = int(nw)
+
  self.stride = stride
- record_t = int(self.nt/stride)
+ record_t = int(self.nt/stride)+1
  # Get starts
- walkers = self._get_starts(same_start)
+ if nw is not None:
+ self.nw = int(nw)
+ walkers = self._get_starts(same_start)
+ else:
+ walkers = np.argwhere(self.im > 0.0)
+ self.nw = walkers.shape[0]
  if self.seed:
  # Generate a seed for each timestep
  np.random.seed(1)
@@ -388,7 +395,7 @@ def plot_msd(self):
  self.data = {}
  fig, ax = plt.subplots(figsize=[6, 6])
  ax.set(aspect=1, xlim=(0, self.nt), ylim=(0, self.nt))
- x = np.arange(0, self.nt, self.stride)[:, np.newaxis]
+ x = np.arange(0, self.nt+1, self.stride)[:, np.newaxis]
  plt.plot(x, self.msd, 'k-', label='msd')
  print('#'*30)
  print('Square Displacement:')
@@ -663,3 +670,14 @@ def run_analytics(self, ws, ts, fname='analytics.csv', num_proc=None):
  header_written = True
  w.writerow(self.data)
  gc.collect()
+
+ def colour_sq_disp(self):
+ starts = self.real_coords[0, :, :]
+ self.im_sq_disp = self.im.copy().astype(int)
+ if self.dim == 3:
+ self.im_sq_disp[starts[:, 0],
+ starts[:, 1],
+ starts[:, 2]] = self.sq_disp[-1, :]
+ else:
+ self.im_sq_disp[starts[:, 0],
+ starts[:, 1]] = self.sq_disp[-1, :]