OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

moved workflows to example directory

Browse Source
This commit is contained in:
James E McClure
2016-11-07 13:04:38 -05:00
parent 74e945c9a3
commit 55df5b0c4d
12 changed files with 839 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
example/CMakeLists.txt
View File

@@ -3,6 +3,10 @@ INSTALL_EXAMPLE( Bubble )
INSTALL_EXAMPLE( ConstrainedBubble )
INSTALL_EXAMPLE( Piston )
INSTALL_EXAMPLE( Sph1896 )
INSTALL_EXAMPLE( drainage )
INSTALL_EXAMPLE( imbibition )
INSTALL_EXAMPLE( relperm )
# Create unit tests for each example
# TEST_EXAMPLE( example exe N_procs )
@@ -12,3 +16,5 @@ INSTALL_EXAMPLE( Sph1896 )

52
example/SegmentingImage/Check_SignDist_and_ID.py Normal file
View File

@@ -0,0 +1,52 @@
#!/usr/bin/env python
import numpy as np
import matplotlib.pylab as plt
#from glob import glob
import sys
# Check if there is a proper command line argument
#if len(sys.argv) !=2:
# sys.stderr.write('Usage: ' + sys.argv[0] + ' Domain.in\n')
# sys.exit()
## end if
# Read 'Domain.in' file
f = open('Domain.in','r') # read-only
lines = f.readlines()
nprocx, nprocy, nprocz = np.fromstring(lines[0].splitlines()[0],dtype=np.int32,sep=' ')
nx, ny, nz = np.fromstring(lines[1].splitlines()[0],dtype=np.int32,sep=' ')
Lx, Ly, Lz = np.fromstring(lines[3].splitlines()[0],dtype=np.float32,sep=' ')
f.close()
lx = 3
ly = 500
lz = 500
# Load the data
SignDist = np.fromfile('SignDist.00000',dtype = np.float64)
SignDist.shape = (lz+2,ly+2,lx+2)
ID = np.fromfile('ID.00000',dtype = np.uint8)
ID.shape = (lz+2,ly+2,lx+2)
# Plot
#plt.figure(1)
#plt.title('SignDist Map')
#plt.pcolormesh(SignDist[:,:,2])
#plt.grid(True)
#plt.axis('equal')
#plt.colorbar()
#
#
#plt.figure(2)
#plt.title('ID.xxxxx')
#plt.pcolormesh(ID[:,:,2],cmap='hot')
#plt.grid(True)
#plt.axis('equal')
#plt.show()

93
example/SegmentingImage/SegmentMicromodelTiff_Group.py Normal file
View File

@@ -0,0 +1,93 @@
#!/usr/bin/env python
import numpy as np
import matplotlib.pylab as plt
from glob import glob
from PIL import Image # import Python Imaging Library (PIL)
import sys
from SegmentMicromodelTiff_utils import *
# Check if there is a proper command line argument
if len(sys.argv) !=2:
sys.stderr.write('Usage: ' + sys.argv[0] + ' <Input parameter file>\n')
print 'Input parameter file: '
print 'line 0: the base name of the experimental images\n',\
'line 1: imin imax\n',\
'line 2: imin_b imax_b\n',\
'line 3: top_layer bottom_layer DEPTH\n',\
'line 4: threshold_nw threshold_s'
sys.exit()
# end if
input_parameter_file = sys.argv[1] # Give the name of the input parameter file
BaseName_output_init_config = '_InitConfig_segmented.raw'# The string will be appended at the end of the
# name of the input image as an output file name
BaseName_output_domain = '_segmented.raw'# The string will be appended at the end of the
# name of the input image as an output file name
image_format = '.tiff' # format of experimental images
# Read other input parameters to get NX, NY and NZ
Para = read_input_parameters(input_parameter_file)
NY = Para['imax'] - Para['imin']
NZ = NY
NX = Para['DEPTH']+Para['top_layer']+Para['bot_layer']
# Load a group of image files with 'base_name' in the names of files
# e.g. 'Micromodel_1.tiff' etc.
file_input_group = glob('*'+Para['base_name']+'*'+image_format) # need to match the image format
# Process all imported experimental images
if not file_input_group:
print 'Error: Input files cannot be found ! '
else:
for ii in range(len(file_input_group)):
file_input_single = file_input_group[ii]
print "Processing image "+file_input_single+" now..."
print "------ Micromodel dimensions (NX, NY, NZ) are (%i, %i, %i) ------" % (NX,NY,NZ)
# Get an array from the input .tiff image
im_array = convert_image_to_array(file_input_single,Para['imin'],Para['imax'])
# Get initial fluid configuration shown in the *.Tiff image
ID = get_initial_config(im_array,Para,NX,NY,NZ)
# Get only the domain (0: solid, 1: fluid)
PorousMedium = get_porous_medium(im_array,Para,NX,NY,NZ)
# Calculate the porosity
POROSITY = 1.0 - 1.0*ID[ID==0].size/ID.size
print "Porosity of the micromodel is: "+str(POROSITY)
# Write the segmeted data to the output
# set the output file names (need to chop off say the '.tiff' part)
output_file_name_domain = file_input_single[:file_input_single.find('.')]+BaseName_output_domain
output_file_name_init_config = file_input_single[:file_input_single.find('.')]+BaseName_output_init_config
# Write out the segmented data (in binary format)
print "The segmented fluid configuration is written to "+output_file_name_init_config+" now..."
print "The segmented porous domain is written to "+output_file_name_domain+" now..."
ID.tofile(output_file_name_init_config)
PorousMedium.tofile(output_file_name_domain)
# NOTE when you want to load the binary file
# Do as follows:
# ID_reload = np.fromfile(file_name,dtype=np.uint8)
#end for
#end if
### In case you want to test the segmentation
#plt.figure(1)
#plt.subplot(1,2,1)
#plt.title('original RGB figure')
#plt.pcolormesh(im_array);
#plt.axis('equal')
#plt.colorbar()
#plt.subplot(1,2,2)
#plt.title('Segmented image')
### This will show the last-processed segmented image
#cmap = plt.cm.get_cmap('hot',3) #Plot 3 discrete colors for NW, W and Solids
#cax=plt.pcolormesh(ID[:,:,NX/2],cmap=cmap,vmin=-0.5,vmax=2.5);
#cbar = plt.colorbar(cax,ticks=[0,1,2])
#plt.axis('equal')
#plt.show()

91
example/SegmentingImage/SegmentMicromodelTiff_utils.py Normal file
View File

@@ -0,0 +1,91 @@
#!/usr/bin/env python
import numpy as np
from PIL import Image # import Python Imaging Library (PIL)
def convert_image_to_array(file_name,imin,imax):
# file_name: the name of the input image (a string)
# imin: (left) boundary that defines the output array
# imax: (right) boundary that defines the output array
# Note: assume a grayscale image (whose R, G, and B values are the same)
# Return an numpy array with RGB values that has the same size as the wanted micromodel
im = Image.open(file_name).convert('RGB')
im_array = np.array(im) # 'im_array' has dimension (height, width, 3)
# Grayscale input image is assumed, thus its R, G, and B values are the same
im_array = im_array[:,:,0] # 2D array with dimension (height, width)
im_array = im_array[:,imin:imax] # Chop the original image to a desirable size
#im_array = np.flipud(im_array[:,imin:imax]) # You might want to flip the Z-axis
return im_array
#end def
def read_input_parameters(input_file_name):
# input_file_name is a string
# The *.in file has the following lines
# line 0: the base name of the experimental images
# line 1: imin imax
# line 2: imin_b imax_b
# line 3: top_layer bottom_layer DEPTH
# line 4: threshold_nw threshold_s
# Note: 'threshold_nw' means: NW phase: RGB values > threshold_nw
# 'threshold_s' means: solid: RGB values < threshold_s
f = open(input_file_name,'r') # read-only
lines = f.readlines()
output_file = {'base_name':lines[0].splitlines()[0]}
line1_array = np.fromstring(lines[1].splitlines()[0],dtype=np.int32,sep=' ')
line2_array = np.fromstring(lines[2].splitlines()[0],dtype=np.int32,sep=' ')
line3_array = np.fromstring(lines[3].splitlines()[0],dtype=np.int32,sep=' ')
line4_array = np.fromstring(lines[4].splitlines()[0],dtype=np.int32,sep=' ')
output_file['imin']=line1_array[0]
output_file['imax']=line1_array[1]
output_file['imin_b']=line2_array[0]
output_file['imax_b']=line2_array[1]
output_file['top_layer']=line3_array[0]
output_file['bot_layer']=line3_array[1]
output_file['DEPTH']=line3_array[2]
output_file['threshold_nw']=line4_array[0]
output_file['threshold_s'] =line4_array[1]
f.close()
return output_file
#end def
def get_initial_config(image_array,Para,NX,NY,NZ):
# This function will give the same fluid configuration shown in the *.Tiff image
# Function input:
# 'image_array' is a numpy array of RGB values, with the same size as micromodel
# 'Para' is a dictionary of input parameters
# Initialize the segmented image 'ID'
# NOTE: 'ID' has the dimension (height, width, DEPTH) or (NZ,NY,NX)
ID = np.zeros((NZ,NY,NX),dtype=np.uint8)
# Map the pixels to the 3D geometry
# Rules: NW phase: RGB values > threshold_nw
# Solid: RGB values < threshold_s
# 1. Identify the non-wetting phase
ID[image_array>Para['threshold_nw'],Para['top_layer']:Para['top_layer']+Para['DEPTH']] = 1
# 2. Identify the wetting phase
ID[np.logical_and(image_array>=Para['threshold_s'],image_array<=Para['threshold_nw']),\
Para['top_layer']:Para['top_layer']+Para['DEPTH']] = 2
# 3. Post boundary retreatment along y-axis. Note z-axis is always the flow direction
ID[:,:Para['imin_b']-Para['imin'],Para['top_layer']:Para['top_layer']+Para['DEPTH']]=0
ID[:,ID.shape[1]-(Para['imax']-Para['imax_b']):,Para['top_layer']:Para['top_layer']+Para['DEPTH']]=0
return ID
#end def
def get_porous_medium(image_array,Para,NX,NY,NZ):
# This function only gives the domain (0->solid,1->fluid) of the *.Tiff image
# Function input:
# 'image_array' is a numpy array of RGB values, with the same size as micromodel
# 'Para' is a dictionary of input parameters
domain = np.zeros((NZ,NY,NX),dtype=np.uint8)
# Map the pixels to the 3D geometry
# Rules: NW phase: RGB values > threshold_nw
# Solid: RGB values < threshold_s
# 1. Identify fluid phase
domain[image_array>=Para['threshold_s'],Para['top_layer']:Para['top_layer']+Para['DEPTH']]=1
# 2. Post boundary retreatment along y-axis. Note z-axis is always the flow direction
domain[:,:Para['imin_b']-Para['imin'],Para['top_layer']:Para['top_layer']+Para['DEPTH']]=0
domain[:,domain.shape[1]-(Para['imax']-Para['imax_b']):,Para['top_layer']:Para['top_layer']+Para['DEPTH']]=0
return domain
#end def

67
example/drainage/eos-MorphDrain.pbs Normal file
View File

@@ -0,0 +1,67 @@
#!/bin/bash
#PBS -A GEO106
#PBS -N MorphDrain
#PBS -j oe
##PBS -l walltime=02:00:00,nodes=216
#PBS -l walltime=01:00:00,nodes=4
##PBS -l gres=widow2%widow3
##PBS -q killable
##PBS -q debug
#cd /tmp/work/$USER
date
cd $PBS_O_WORKDIR
#LBPM_WIA_INSTALL_DIR=/lustre/atlas/proj-shared/geo106/build-eos-LBPM-WIA
LBPM_WIA_INSTALL_DIR=$MEMBERWORK/geo106/eos-LBPM-WIA
#echo "PBS_O_WORKDIR: `echo $PBS_O_WORKDIR`"
source $MODULESHOME/init/bash
module swap PrgEnv-intel PrgEnv-gnu
#module swap cray-mpich2 cray-mpich2/5.6.3
#module load cudatoolkit
export LD_LIBRARY_PATH=${CRAY_LD_LIBRARY_PATH}:${LD_LIBRARY_PATH}
#export MPICH_RDMA_ENABLED_CUDA=1
#aprun -n 27 -N 1 /lustre/atlas/scratch/mcclurej/geo106/LBPM-WIA/bin/lb2_Color_wia_mpi
#LIST=$(ls|grep sw)
#for DIR in $LIST; do
# cd $DIR
# sat=$(cat Color.in |head -3 | tail -1)
# aprun -n 64 $LBPM_WIA_INSTALL_DIR/bin/lbpm_segmented_decomp
# aprun -n 64 $LBPM_WIA_INSTALL_DIR/bin/lbpm_segmented_pp
# aprun -n 144 $LBPM_WIA_INSTALL_DIR/bin/lbpm_random_pp 1 0
# aprun -n 144 $LBPM_WIA_INSTALL_DIR/bin/lbpm_random_pp 1 1
# aprun -n 64 $LBPM_WIA_INSTALL_DIR/bin/lbpm_segmented_decomp 1 2
# aprun -n 64 $LBPM_WIA_INSTALL_DIR/bin/lbpm_segmented_pp
# aprun -n 64 $LBPM_WIA_INSTALL_DIR/bin/lbpm_random_pp 0 1
#mkdir -p MEDIA
#cp ID* MEDIA/
cp MEDIA/ID* ./
LABEL="lrc32_sandpack"
# Run morphological drainage and set up input files for media
RADIUS="5.4 5.35 5.2 5 4.8 4.6 4.4 4.2 4 3.67 3.33 3 2.67 2.33 2 1.5"
echo "radius sw" > morphdrain.csv
for r in $RADIUS; do
aprun -n 64 $LBPM_WIA_INSTALL_DIR/bin/lbpm_morphdrain_pp $r > morph.log
sat=$(grep sat morph.log | sed 's/Final saturation=//g')
tag=$(wc -l morphdrain.csv | cut -c1-2)
echo "$r $sat" >> morphdrain.csv
DIR=$LABEL"_drain_"$tag
mkdir -p $DIR
cp ID.* $DIR
cp SignDist.* $DIR
cp Domain.in $DIR
done
exit;

74
example/drainage/setup-drain.py Executable file
View File

@@ -0,0 +1,74 @@
import os
import sys
import csv
name=sys.argv[1]
#numpts=int(sys.argv[2])
process="drain"
cwd=os.getcwd()
print("Initializing morphological drainage cases")
print("Base name is "+name)
# Parameters for color LBM
tau=0.7
alpha=0.005
beta=0.95
phisol=-1.0
saturation=0.0
Fx=0.0
Fy=0.0
Fz=0.0
Restart=0
pBC=1
din=1.001
dout=0.999
maxtime=100005
interval=50000
tolerance=1e-5
viscosity=(tau-0.5)/3
ift=5.796*alpha
radius=[]
sw=[]
with open("morphdrain.csv","r") as f:
for line in f:
reader=csv.reader(f,delimiter=' ')
for row in reader:
radius.append(float(row[0]))
sw.append(float(row[1]))
numpts=len(radius)
print("Number of cases "+str(numpts))
for pt in range(0,numpts):
#compute the pressure difference
tag=pt+1
dp=2*ift/radius[pt]
din=1.0+0.5*dp
dout=1.0-0.5*dp
dirname=str(name)+"_"+str(process)+"_"+str(tag)
print("Creating " + dirname)
if not os.path.exists(dirname):
os.mkdir(dirname)
os.chdir(dirname)
ParamFile = open("Color.in","w")
ParamFile.write("%f\n" % tau)
ParamFile.write("%f " % alpha)
ParamFile.write("%f " % beta)
ParamFile.write("%f\n" % phisol)
ParamFile.write("%f\n" % saturation)
ParamFile.write("%f " % Fx)
ParamFile.write("%f " % Fy)
ParamFile.write("%f\n" % Fz)
ParamFile.write("%i " % Restart)
ParamFile.write("%i " % pBC)
ParamFile.write("%f " % din)
ParamFile.write("%f\n" % dout)
ParamFile.write("%i " % maxtime)
ParamFile.write("%i " % interval)
ParamFile.write("%f\n" % tolerance)
ParamFile.close()
os.chdir(cwd)

20
example/imbibition/InitializeImbibition.sh Executable file
View File

@@ -0,0 +1,20 @@
#!/bin/bash
# Run morphological analysis and set up input files for media
LABEL="sph1964"
#aprun -n 64 $LBPM_WIA_INSTALL_DIR/bin/lbpm_morph_pp $r > morph.log
echo "Quartile Radius" > poresize.quartiles
grep -A4 Quartile morph.log | tail -4 >> poresize.quartiles
python ~/LBPM-WIA/workflows/imbibition/setup-imbibition.py $LABEL
FILES=$(ls | grep "drain_")
for src in $FILES; do
dest=$(echo $src | sed 's/drain_/imb_/g')
echo "Initializing $dest from primary drainage"
cp -r $src $dest
cp Color.in.imbibition $dest/Color.in
done
exit;

77
example/imbibition/setup-imbibition.py Executable file
View File

@@ -0,0 +1,77 @@
import os
import sys
import csv
import glob
name=sys.argv[1]
#numpts=int(sys.argv[2])
process="drain"
cwd=os.getcwd()
print("Initializing imbibition sequence")
print("Drainage should initialize imbibition!")
print("Base name is "+name)
# Parameters for color LBM
tau=0.8
alpha=0.005
beta=0.95
phisol=-1.0
saturation=0.0
Fx=0.0
Fy=0.0
Fz=0.0
Restart=1
pBC=3
din=1.001
dout=0.999
maxtime=2000005
interval=50000
tolerance=1e-5
viscosity=(tau-0.5)/3
ift=5.796*alpha
radius=[]
#porsize.quartiles provides quartiles for the pore size
with open("poresize.quartiles","r") as f:
for line in f:
reader=csv.reader(f,delimiter=' ')
for row in reader:
radius.append(float(row[1]))
#sw.append(float(row[1]))
#compute the pressure difference
# Use Q1 as the initial pressure
din=2*ift/radius[0]
# Use Q4 (maximum pore size) to set final pressure differnece
dout=2*ift/(radius[3]+1.0)
# iterator for the filenames matching drain in the current directory
#Source=glob.iglob('*drain_*')
#dirname=str(name)+"_"+str(process)+"_"+str(pt)
#print("Creating " + dirname)
#if not os.path.exists(dirname):
# os.mkdir(dirname)
#os.chdir(dirname)
ParamFile = open("Color.in.imbibition","w")
ParamFile.write("%f\n" % tau)
ParamFile.write("%f " % alpha)
ParamFile.write("%f " % beta)
ParamFile.write("%f\n" % phisol)
ParamFile.write("%f\n" % saturation)
ParamFile.write("%f " % Fx)
ParamFile.write("%f " % Fy)
ParamFile.write("%f\n" % Fz)
ParamFile.write("%i " % Restart)
ParamFile.write("%i " % pBC)
ParamFile.write("%f " % din)
ParamFile.write("%f\n" % dout)
ParamFile.write("%i " % maxtime)
ParamFile.write("%i " % interval)
ParamFile.write("%f\n" % tolerance)
ParamFile.close()
#os.chdir(c)

270
example/relperm/dist_func_utils.py Normal file
View File

@@ -0,0 +1,270 @@
#!/usr/bin/env python
import numpy as np
from scipy import ndimage
from scipy import spatial
import scipy.ndimage.filters as filters
import scipy.ndimage.morphology as morphology
import scipy.stats as stats
from skimage.morphology import medial_axis,skeletonize_3d
def detect_intersection_point_with_cluster_filter(input_data,lx,ly,lz=0):
# Input is a signed distance data file such as 'SignDist.xxxxx'
if lz > 0: # i.e. a 3D input
# Calculate the medial axis of the segmented image
dist = input_data.copy()
dist.shape = (lz,ly,lx)
skel = skeletonize_3d(dist>0.0)
dist_on_skel = skel*dist
[z_skel,y_skel,x_skel]=np.where(skel)
# Building two search trees is potentially good for 3D image
grid = np.indices(dist.shape)
grid_points = zip(grid[0].ravel(),grid[1].ravel(),grid[2].ravel())
tree_grid = spatial.cKDTree(grid_points)
points_for_search = zip(z_skel,y_skel,x_skel)
tree_points_for_search = spatial.cKDTree(points_for_search)
neighbor_all = tree_points_for_search.query_ball_tree(tree_grid,np.sqrt(3.0))
idx_glb_table = np.ravel_multi_index([z_skel,y_skel,x_skel],skel.shape)
idx_glb_candidate = np.empty((0,),dtype=int)
for k,idx_glb_neighbor in enumerate(neighbor_all):
if k%4000==0:
print 'Search for intersection points: '+str(k)+'/'+str(len(neighbor_all)-1)
#end if
mask = np.in1d(idx_glb_neighbor,idx_glb_table,assume_unique=True)
idx_glb_candidate = np.hstack((idx_glb_candidate,mask.sum()))
#end for
# Statistics: number of neighbors for each voxel on the medial axis
connectivity_stats = stats.itemfreq(idx_glb_candidate)
# 'connectivity_stats' has the following format:
# number_of_neighbors_for_each_voxel :: corresponding_number_of_voxels
# array([[ 1, 41],
# [ 2, 143],
# [ 3, 185],
# [ 4, 5]])
#
# If a voxel is indentified as an intersection point, the number of its
# neighboring voxels should be greater than 'benchmark'
benchmark = connectivity_stats[np.argmax(connectivity_stats[:,1]),0]
# Update the medial axis and 'dist_on_skel'
skel[:] = False
skel[np.unravel_index(idx_glb_table[idx_glb_candidate > benchmark],skel.shape)] = True
else: # i.e. 2D input
# Calculate the medial axis of the segmented image
dist = input_data.copy()
dist.shape = (ly,lx)
skel = medial_axis(dist>0.0)
dist_on_skel = skel*dist
[y_skel,x_skel]=np.where(skel)
# Building two search trees is potentially good for 3D image
grid = np.indices(dist.shape)
grid_points = zip(grid[0].ravel(),grid[1].ravel())
tree_grid = spatial.cKDTree(grid_points)
points_for_search = zip(y_skel,x_skel)
tree_points_for_search = spatial.cKDTree(points_for_search)
neighbor_all = tree_points_for_search.query_ball_tree(tree_grid,np.sqrt(2.0))
idx_glb_table = np.ravel_multi_index([y_skel,x_skel],skel.shape)
idx_glb_candidate = np.empty((0,),dtype=int)
for k,idx_glb_neighbor in enumerate(neighbor_all):
if k%4000==0:
print 'Search for intersection points: '+str(k)+'/'+str(len(neighbor_all)-1)
#end if
mask = np.in1d(idx_glb_neighbor,idx_glb_table,assume_unique=True)
idx_glb_candidate = np.hstack((idx_glb_candidate,mask.sum()))
#end for
# Statistics: number of neighbors for each voxel on the medial axis
connectivity_stats = stats.itemfreq(idx_glb_candidate)
# 'connectivity_stats' has the following format:
# number_of_neighbors_for_each_voxel :: corresponding_number_of_voxels
# array([[ 1, 41],
# [ 2, 143],
# [ 3, 185],
# [ 4, 5]])
#
# If a voxel is indentified as an intersection point, the number of its
# neighboring voxels should be greater than 'benchmark'
benchmark = connectivity_stats[np.argmax(connectivity_stats[:,1]),0]
# Update the medial axis and 'dist_on_skel'
skel[:] = False
skel[np.unravel_index(idx_glb_table[idx_glb_candidate > benchmark],skel.shape)] = True
#end if
return (_Filter_cluster_close_points(skel*dist),connectivity_stats)
#end def
def _Filter_cluster_close_points(arr):
# 'arr' is a 2D/3D signed distance data file
# Return an 'arr' with nearest neighboring points clustered
if arr.ndim == 2:
[y_idx,x_idx] = np.where(arr>0.0)
idx_glb = np.ravel_multi_index([y_idx,x_idx],arr.shape)
grid = np.indices(arr.shape)
dist = arr.ravel()[idx_glb]
candidate = np.empty((0,),dtype=int)
# TODO: use search tree to do this !
for k in np.arange(idx_glb.size):
if k%200==0:
print 'Clustering close points: '+str(k)+'/'+str(idx_glb.size-1)
#end if
mask_circle = (grid[0]-y_idx[k])**2 + (grid[1]-x_idx[k])**2<=dist[k]**2
idx_glb_circle =np.ravel_multi_index(np.where(mask_circle),mask_circle.shape)
mask = np.in1d(idx_glb,idx_glb_circle,assume_unique=True)
candidate = np.hstack((candidate,idx_glb[mask][np.argmax(dist[mask])]))
#end for
candidate = np.unique(candidate)
mask = np.in1d(idx_glb,candidate,assume_unique=True)
output_glb_idx = idx_glb[mask]
output = np.zeros_like(arr.ravel())
output[output_glb_idx] = arr.ravel()[output_glb_idx]
output.shape = arr.shape
else:
[z_idx,y_idx,x_idx] = np.where(arr>0.0)
idx_glb = np.ravel_multi_index([z_idx,y_idx,x_idx],arr.shape)
grid = np.indices(arr.shape)
dist = arr.ravel()[idx_glb]
candidate = np.empty((0,),dtype=int)
# TODO: use search tree to do this !
for k in np.arange(idx_glb.size):
if k%200==0:
print 'Clustering close points: '+str(k)+'/'+str(idx_glb.size-1)
#end if
mask_circle = (grid[0]-z_idx[k])**2 + (grid[1]-y_idx[k])**2 + (grid[2]-x_idx[k])**2<=dist[k]**2
idx_glb_circle =np.ravel_multi_index(np.where(mask_circle),mask_circle.shape)
mask = np.in1d(idx_glb,idx_glb_circle,assume_unique=True)
candidate = np.hstack((candidate,idx_glb[mask][np.argmax(dist[mask])]))
#end for
candidate = np.unique(candidate)
mask = np.in1d(idx_glb,candidate,assume_unique=True)
output_glb_idx = idx_glb[mask]
output = np.zeros_like(arr.ravel())
output[output_glb_idx] = arr.ravel()[output_glb_idx]
output.shape = arr.shape
#end if
return output
#end def
def get_Dist(f):
return ndimage.distance_transform_edt(f)
#end def
def detect_local_maxima(arr,medial_axis_arr,patch_size=3):
local_max = _find_local_maxima(arr,patch_size)
background_value = 1e5 # for denoising in finding local minima
arr2 = arr.copy()
arr2[np.logical_not(medial_axis_arr)]=background_value
ocal_min = _find_local_minima(arr2,background_val=background_value)
local_min = np.logical_and(medial_axis_arr,local_min) #Correct min_indices with Medial axis
local_max_reduced = np.logical_xor(local_max,local_min)
local_max = np.logical_and(local_max,local_max_reduced)
local_max = np.logical_and(medial_axis_arr,local_max)#Correct max_indices with Medial axis
return local_max
#end def
def detect_local_maxima_with_cluster_filter(arr,medial_axis_arr,patch_size=3):
# apply the cluster filetering only once
local_max = detect_local_maxima(arr,medial_axis_arr,patch_size=patch_size)
arr2 = arr.copy()
arr2[np.logical_not(local_max)]=0.0
return _Filter_cluster_close_points(arr2)
#end def
def detect_local_maxima_with_cluster_filter_loop(arr,medial_axis_arr,patch_start=3,patch_stop=9):
# Implement multiple search for local maxima (with cluster filtering)
arr1 = detect_local_maxima_with_cluster_filter(arr,medial_axis_arr,patch_size=patch_start)
[y_idx_1,x_idx_1] = np.where(arr1>0.0)
arr2 = detect_local_maxima_with_cluster_filter(arr1,medial_axis_arr,patch_size=patch_start+2)
[y_idx_2,x_idx_2] = np.where(arr2>0.0)
if (y_idx_1.size>y_idx_2.size):
counter = 2 # Record how many times is the filtering applied
patch = np.arange((patch_start+2)+2,patch_stop+2,2)
for k in patch:
y_idx_1 = y_idx_2.copy()
x_idx_1 = x_idx_2.copy()
arr1 = detect_local_maxima_with_cluster_filter(arr2,medial_axis_arr,patch_size=k)
[y_idx_2,x_idx_2] = np.where(arr1>0.0)
arr2 = arr1.copy()
counter = counter + 1
if not (y_idx_1.size > y_idx_2.size):
print 'Note: Local maxima are already found before the patch_stop is reached.'
print ' All local maxima are found at the patch size of '+str(k)+' !'
break;
#end if
#end for
else:
print 'Note: All local maxima are found at the patch size of '+str(patch_start)+' !'
return arr2
#end if
return arr2
#end def
#def Filter_cluster(arr):
# return _Filter_cluster_close_points(arr)
##end def
def _find_local_maxima(arr,patch_size):
# http://stackoverflow.com/questions/3684484/peak-detection-in-a-2d-array/3689710#3689710
"""
Takes an array and detects the troughs using the local maximum filter.
Returns a boolean mask of the troughs (i.e. 1 when
the pixel's value is the neighborhood maximum, 0 otherwise)
"""
# define an connected neighborhood
# http://www.scipy.org/doc/api_docs/SciPy.ndimage.morphology.html#generate_binary_structure
#neighborhood = morphology.generate_binary_structure(len(arr.shape),2)
#neighborhood = np.ones((patch_size,patch_size),dtype=bool)
# apply the local maximum filter; all locations of maximum value
# in their neighborhood are set to 1
local_max = (filters.maximum_filter(arr, size=patch_size,mode='constant',cval=0)==arr)
# local_min is a mask that contains the peaks we are
# looking for, but also the background.
# In order to isolate the peaks we must remove the background from the mask.
#
# we create the mask of the background
background = (arr==0)
#
# a little technicality: we must erode the background in order to
# successfully subtract it from local_min, otherwise a line will
# appear along the background border (artifact of the local minimum filter)
# http://www.scipy.org/doc/api_docs/SciPy.ndimage.morphology.html#binary_erosion
if arr.ndim == 3:
neighborhood = np.ones((patch_size,patch_size,patch_size),dtype=bool)
elif arr.ndim==2:
neighborhood = np.ones((patch_size,patch_size),dtype=bool)
else:
neighborhood = np.ones((patch_size,),dtype=bool)
#end if
eroded_background = morphology.binary_erosion(
background, structure=neighborhood, border_value=1)
#
# we obtain the final mask, containing only peaks,
# by removing the background from the local_min mask
detected_maxima = local_max - eroded_background
return detected_maxima
def _find_local_minima(arr,patch_size=3,background_val=1e5):
local_min = (filters.minimum_filter(arr, size=patch_size,mode='constant',cval=background_val)==arr)
background = (arr==background_val)
if arr.ndim == 3:
neighborhood = np.ones((patch_size,patch_size,patch_size),dtype=bool)
elif arr.ndim==2:
neighborhood = np.ones((patch_size,patch_size),dtype=bool)
else:
neighborhood = np.ones((patch_size,),dtype=bool)
#end if
eroded_background = morphology.binary_erosion(
background, structure=neighborhood, border_value=1)
detected_minima = local_min - eroded_background
return detected_minima
#end def

56
example/relperm/morph_getpores.py Normal file
View File

@@ -0,0 +1,56 @@
#!/usr/bin/env python
import sys
import numpy as np
from dist_func_utils import *
from glob import glob
# Check if there is a proper command line argument
if len(sys.argv) !=2:
sys.stderr.write('Usage: ' + sys.argv[0] + ' <Domain.in>\n')
sys.exit()
# end if
# Read 'Domain.in' to obtain the size of 'SignDist.xxxxx'
f = open(sys.argv[1],'r')
lines = f.readlines()
nx,ny,nz = np.fromstring(lines[1].splitlines()[0],dtype = np.int32,sep=' ')
f.close()
#nx = ny = nz = 128
nx+=2 # consider the halo layer
ny+=2
nz+=2
base_name = 'SignDist.'
file_input_group = glob(base_name+'*')
# Prepare output file: 'pores_xxxxx.csv'
output_data_name = 'pores_'
output_data_format = '.csv'
# Process all imported experimental images
if not file_input_group:
print 'Error: Input files cannot be found ! '
else:
for ii in range(len(file_input_group)):
file_input_single = file_input_group[ii]
file_input_single_idx = file_input_single[file_input_single.find(base_name)+len(base_name):]
print '--- Get pore size information for '+file_input_single+' now ---'
dist = np.fromfile(file_input_single,dtype = np.float64)
dist_on_skel,connect_stats = detect_intersection_point_with_cluster_filter(dist,nx,ny,nz)
[z_skel,y_skel,x_skel] = np.where(dist_on_skel>0.0)
output_data = np.column_stack((x_skel,y_skel,z_skel,dist_on_skel[dist_on_skel>0.0]))
print '--- Save pore size csv file ---'
np.savetxt(output_data_name+file_input_single_idx+output_data_format,output_data,delimiter=' ')
#end for
#end if

32
example/relperm/morph_reducepores.py Normal file
View File

@@ -0,0 +1,32 @@
#!/usr/bin/env python
import numpy as np
import sys
# Check if there is a proper command line argument
if len(sys.argv) !=3:
sys.stderr.write('Usage: ' + sys.argv[0] + ' pores.csv Sw_list\n')
sys.exit()
# end if
# Read 'pores.csv' and a list of Sw
pores = np.genfromtxt(sys.argv[1])
Sw = np.genfromtxt(sys.argv[2])
#NOTE: 'pores.csv' has a layout of : [x,y,z,pore_radius]
pores = pores[:,-1]
# Calculate the percentile according to the Sw list
if Sw.max()<=1.0 and Sw.min()>=0.0:
radii_init_imbib = np.percentile(pores,100.0*Sw.ravel())
else:
print 'Error: the list of Sw should have values 0.0 - 1.0 !'
sys.exit()
#end if
radii_init_imbib.shape = (radii_init_imbib.size,1)
# Write out initial guess for imbibition
output_name = 'radius.csv'
print '------ Initial radii for morphological imbibition is written to the file: '+output_name+' ------'
np.savetxt(output_name,radii_init_imbib)