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

add morphological drainage and imbibition

Browse Source
This commit is contained in:
zherexli
2016-11-01 15:57:27 -04:00
parent 73ce928dd2
commit 13a533a748
2 changed files with 516 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

119
workflows/Morphological_Analysis/generate_morphological_configs.py Normal file
View File

@@ -0,0 +1,119 @@
#!/usr/bin/env python
import numpy as np
import matplotlib.pylab as plt
from glob import glob
import sys
#from dist_func_utils import *
import scipy.stats as stats
import scipy.ndimage.morphology as morphology
import scipy.ndimage.measurements as measurements
from morphological_analysis_utils import *
# Check if there is a proper command line argument
#if len(sys.argv) !=2:
# sys.stderr.write('Usage: ' + sys.argv[0] + ' <Base name of the segmented image data>\n')
# sys.exit()
# end if
base_name = 'MicroModel' # Give the name of the input parameter file
#base_name = 'SignDist' # Give the name of the input parameter file
#base_name = sys.argv[1] # Give the name of the input parameter file
SegImage_format = '.raw'
# Load a group of segmented image data files with 'base_name' in the names of files
# e.g. 'Micromodel_1_segmented.raw' etc.
file_input_group = glob('*'+base_name+'*'+SegImage_format) # need to match the image data format
# Dimensions for segmented image
lz = 500
ly = 500
lx = 12
#lz = 200
#ly = 200
#lx = 200
R_critical = 4.5
# 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 "Analyzing segmented image data "+file_input_single+" now..."
# The imported data is a double-precision signed distance map
image_raw = np.fromfile(file_input_single,dtype=np.uint8)
#image_raw = np.fromfile(file_input_single,dtype=np.float64)
image_raw.shape=(lz,ly,lx)
edt_dist = morphology.distance_transform_edt(image_raw)
#edt_dist = edt_dist[:,:,5]
# wrap the medium with one layer of solid phase
#image_raw = np.lib.pad(image_raw,((0,0),(1,1),(1,1)),'constant',constant_values=0.0)
drain_config,Sw_drain = generate_morph_drain_config(edt_dist,R_critical)
imbib_config,Sw_imbib = generate_morph_imbib_config(edt_dist,R_critical)
#drain_config,Sw_drain = generate_morph_drain_config(image_raw,R_critical)
#imbib_config,Sw_imbib = generate_morph_imbib_config(image_raw,R_critical)
#Sw_drain = generate_morph_drain_curv_3D(image_raw,R_critical)
#imbib_config,Sw_imbib = generate_morph_imbib_config2_3D(image_raw,R_critical)
#Sw_imbib = generate_morph_imbib_curv2_3D(image_raw,R_critical)
print 'Rcrit: '+str(R_critical)+', Drain_Sw = '+str(Sw_drain)
print 'Rcrit: '+str(R_critical)+', Imbib_Sw = '+str(Sw_imbib)
# Save the configuration files
# drain_config.tofile('Drain_config_Rcrit_'+str(R_critical)+'.raw')
# imbib_config.tofile('Imbib_config_Rcrit_'+str(R_critical)+'.raw')
# Load the saved data
# drain_config = np.fromfile('Drain_config_Rcrit_'+str(R_critical)+'.raw',dtype=np.uint8)
# imbib_config = np.fromfile('Imbib_config_Rcrit_'+str(R_critical)+'.raw',dtype=np.uint8)
# drain_config.shape = (lz,ly,lx)
# imbib_config.shape = (lz,ly,lx)
plt.figure(1)
plt.subplot(1,2,1)
plt.title('Drainage: Rcrit='+str(R_critical)+' Sw='+str(Sw_drain))
plt.pcolormesh(drain_config[:,:,lx/2],cmap='hot')
#plt.pcolormesh(drain_config[20,:,:],cmap='hot')
plt.axis('equal')
plt.grid(True)
plt.subplot(1,2,2)
plt.title('Imbibition: Rcrit='+str(R_critical)+' Sw='+str(Sw_imbib))
plt.pcolormesh(imbib_config[:,:,lx/2],cmap='hot')
#plt.pcolormesh(imbib_config[20,:,:],cmap='hot')
plt.axis('equal')
plt.grid(True)
# plt.figure(1)
# plt.subplot(1,2,1)
# plt.title('Drainage: Rcrit='+str(R_critical)+' Sw='+str(Sw_drain))
# plt.pcolormesh(drain_config,cmap='hot')
# plt.axis('equal')
# plt.grid(True)
#
# plt.subplot(1,2,2)
# plt.title('Imbibition: Rcrit='+str(R_critical)+' Sw='+str(Sw_imbib))
# plt.pcolormesh(imbib_config,cmap='hot')
# plt.axis('equal')
# plt.grid(True)
plt.figure(2)
plt.plot(Sw_drain,1.0/R_critical,'ro',markersize=6,label='Drainage')
plt.plot(Sw_imbib,1.0/R_critical,'b^',markersize=6,label='Imbibition')
plt.legend(loc='best')
plt.grid(True)
plt.xlim([0,1.0])
plt.ylim([0,2.0])
plt.show()
#end if
#end for

397
workflows/Morphological_Analysis/morphological_analysis_utils.py Normal file
View File

@@ -0,0 +1,397 @@
#!/usr/bin/env python
import numpy as np
import scipy.stats as stats
import scipy.ndimage.morphology as morphology
import scipy.ndimage.measurements as measurements
#from dist_func_utils import *
def generate_morph_imbib_config(seg_image_input,R_critical):
# This funciton works well at high-Pc-low-Sw
# There is no retarded Sw when the curvature is high
# The method of this functions:
# 1. Perform erosion with the radius of R_critical
# 2. Perform dilation with the radius of R_critical
if seg_image_input.ndim ==2 :
seg_image = seg_image_input>0.0
pore_vol = 1.0*seg_image.sum()
radius = R_critical
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion with radius of R_critical
seg_image_ero = morphology.binary_erosion(seg_image,structure=circle,border_value=1)
# NOTE: Be careful with the 'border_value' of the erosion - should be 'border_value=1'
# Step 2: Perform dilation with radius of R_critical
seg_image_dil = morphology.binary_dilation(seg_image_ero,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of the array after dilation is 'bool'
seg_image_dil = seg_image_dil.astype(np.uint8)
seg_image_dil[np.logical_not(seg_image_dil)]=2
seg_image_dil[np.logical_not(seg_image)] = 0
Sw = (seg_image_dil==2).sum()/pore_vol
else: # 3D
seg_image = seg_image_input>0.0
pore_vol = 1.0*seg_image.sum()
radius = R_critical
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 + (grid[2]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion with radius of R_critical
seg_image_ero = morphology.binary_erosion(seg_image,structure=circle,border_value=1)
# NOTE: Be careful with the 'border_value' of the erosion - should be 'border_value=1'
# Step 2: Perform dilation with radius of R_critical
seg_image_dil = morphology.binary_dilation(seg_image_ero,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of the array after dilation is 'bool'
seg_image_dil = seg_image_dil.astype(np.uint8)
seg_image_dil[np.logical_not(seg_image_dil)]=2
seg_image_dil[np.logical_not(seg_image)] = 0
Sw = (seg_image_dil==2).sum()/pore_vol
#end if
return (seg_image_dil,Sw)
#end def
def generate_morph_imbib_curv(seg_image_input,R_critical):
# This funciton works well at high-Pc-low-Sw
# There is no retarded Sw when the curvature is high
# The method of this functions:
# 1. Perform erosion with the radius of R_critical
# 2. Perform dilation with the radius of R_critical
# *********************************************************************
# Input: seg_image: a well shaped segmented image with size (lz,ly,lx)
# seg_image has values as : NW phase -> 1
# W phase -> 2
# solid phase -> 0
# *********************************************************************
if seg_image_input.ndim == 2:
pore_vol = 1.0*(seg_image_input>0.0).sum()
radius = R_critical
print 'Morphological Imbibition: processing critical radius: '+str(radius)+' now......'
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion with radius of R_critical
seg_image_ero = morphology.binary_erosion(seg_image_input>0.0,structure=circle,border_value=1)
# NOTE: Be careful with the 'border_value' of the erosion - should be 'border_value=1'
# Step 2: Perform dilation with radius of R_critical
seg_image_ero_dil = morphology.binary_dilation(seg_image_ero,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of the array after dilation is 'bool'
seg_image_ero_dil[seg_image_input<=0.0]=False
Sw = 1.0 - seg_image_ero_dil.sum()/pore_vol
else:
pore_vol = 1.0*(seg_image_input>0.0).sum()
radius = R_critical
print 'Morphological Imbibition: processing critical radius: '+str(radius)+' now......'
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 + (grid[2]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion with radius of R_critical
seg_image_ero = morphology.binary_erosion(seg_image_input>0.0,structure=circle,border_value=1)
# NOTE: Be careful with the 'border_value' of the erosion - should be 'border_value=1'
# Step 2: Perform dilation with radius of R_critical
seg_image_ero_dil = morphology.binary_dilation(seg_image_ero,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of the array after dilation is 'bool'
seg_image_ero_dil[seg_image_input<=0.0]=False
Sw = 1.0 - seg_image_ero_dil.sum()/pore_vol
#end if
return Sw
#end def
def generate_morph_drain_config(seg_image_input,R_critical):
# The method for this function follows Hilper & Miller AWR(2001)
# 1. Perform erosion for the pore space with radius of R_critical
# 2. Label the eroded pore space, and leave only the pore space that is still
# connected with the non-wetting phase reservoir
# 3. Perform the dilation for the labelled pore space with radius of R_critical
# **************************************************************************
# Input: seg_image: a well shaped segmented image with size (lz,ly,lx)
# seg_image has values as : NW phase -> 1
# W phase -> 2
# solid phase -> 0
# **************************************************************************
if seg_image_input.ndim ==2:
seg_image = seg_image_input>0.0
pore_vol = 1.0*seg_image.sum()
radius = R_critical
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion on the pore space
# NOTE: the dtype of 'seg_im_ero' is 'bool'
seg_im_ero = morphology.binary_erosion(seg_image,structure=circle,border_value=1)
# NOTE: 'border_value' for erosion should be 'True'
# Step 2: Label the eroded pore space
# NOTE: Assume the NW phase reservoir is at the first layer of the domain
# i.e. at seg_image[0,:] - adjust it if this does not suit your need
# For erosion, assume that diagonals are not considered
# For erosion, assume that diagonals are not considered
seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(2,1))
#seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(2,2))
# NOTE: Here I assume the inlet is at the first layer of the array's axis=2 (i.e. domain[0,:,:])\
# You can always change to any other layers as the inlet for this drainage.
label_check = seg_im_ero_label_temp[0,seg_im_ero_label_temp[0,:]!=0]
label_check = np.unique(label_check)
# NOTE the following lines are only for your to check things
# ******************** For check *******************************#
# It assign the labelled array as: NW -> 1, W -> 2, Solid -> 0
#seg_im_ero_label_show = seg_im_ero_label.copy()
#seg_im_ero_label_show[seg_im_ero_label_show !=1] = 2
#seg_im_ero_label_show[np.logical_not(seg_image_2d)]=0
# ******************** End: for check **************************#
seg_im_ero_label = np.zeros_like(seg_im_ero_label_temp,dtype=bool)
for labels in label_check:
seg_im_ero_label = np.logical_or(seg_im_ero_label,seg_im_ero_label_temp==labels)
seg_im_ero_label = seg_im_ero_label.astype(np.uint8)
# Step 3: perform dilation on the labelled pore space
seg_im_ero_label_dil = morphology.binary_dilation(seg_im_ero_label,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of 'seg_im_ero_label_dil' is 'bool'
seg_im_ero_label_dil = seg_im_ero_label_dil.astype(np.uint8)
seg_im_ero_label_dil[np.logical_not(seg_im_ero_label_dil)]=2
seg_im_ero_label_dil[np.logical_not(seg_image)]=0
Sw = (seg_im_ero_label_dil==2).sum()/pore_vol
else: # 3D porous medium
seg_image = seg_image_input>0.0
pore_vol = 1.0*seg_image.sum()
radius = R_critical
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 + (grid[2]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion on the pore space
# NOTE: the dtype of 'seg_im_ero' is 'bool'
seg_im_ero = morphology.binary_erosion(seg_image,structure=circle,border_value=1)
# NOTE: 'border_value' for erosion should be 'True'
# Step 2: Label the eroded pore space
# NOTE: Assume the NW phase reservoir is at the first layer of the domain
# i.e. at seg_image[0,:] - adjust it if this does not suit your need
# For erosion, assume that diagonals are not considered
seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(3,1))
#seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(3,3))
# NOTE: Here I assume the inlet is at the first layer of the array's axis=2 (i.e. domain[0,:,:])\
# You can always change to any other layers as the inlet for this drainage.
label_check = seg_im_ero_label_temp[0,seg_im_ero_label_temp[0,:]!=0]
label_check = np.unique(label_check)
# NOTE the following lines are only for your to check things
# ******************** For check *******************************#
# It assign the labelled array as: NW -> 1, W -> 2, Solid -> 0
#seg_im_ero_label_show = seg_im_ero_label.copy()
#seg_im_ero_label_show[seg_im_ero_label_show !=1] = 2
#seg_im_ero_label_show[np.logical_not(seg_image_2d)]=0
# ******************** End: for check **************************#
seg_im_ero_label = np.zeros_like(seg_im_ero_label_temp,dtype=bool)
for labels in label_check:
seg_im_ero_label = np.logical_or(seg_im_ero_label,seg_im_ero_label_temp==labels)
seg_im_ero_label = seg_im_ero_label.astype(np.uint8)
# Step 3: perform dilation on the labelled pore space
seg_im_ero_label_dil = morphology.binary_dilation(seg_im_ero_label,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of 'seg_im_ero_label_dil' is 'bool'
seg_im_ero_label_dil = seg_im_ero_label_dil.astype(np.uint8)
seg_im_ero_label_dil[np.logical_not(seg_im_ero_label_dil)]=2
seg_im_ero_label_dil[np.logical_not(seg_image)]=0
Sw = (seg_im_ero_label_dil==2).sum()/pore_vol
#end if
return (seg_im_ero_label_dil,Sw)
#end def
def generate_morph_drain_curv(seg_image_input,R_critical):
# The method for this function follows Hilper & Miller AWR(2001)
# 1. Perform erosion for the pore space with radius of R_critical
# 2. Label the eroded pore space, and leave only the pore space that is still
# connected with the non-wetting phase reservoir
# 3. Perform the dilation for the labelled pore space with radius of R_critical
# ****************************************************************************
# Currently I am provided with a 3D SignDist image which has positive values
# in the pore space and 0.0 in the solid phase.
# ****************************************************************************
if seg_image_input.ndim == 2:
pore_vol = 1.0*(seg_image_input>0.0).sum()
radius = R_critical
print 'Morphological Drainage: processing critical radius: '+str(radius)+' now......'
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion on the pore space
# NOTE: the dtype of 'seg_im_ero' is 'bool'
seg_im_ero = morphology.binary_erosion(seg_image_input>0.0,structure=circle,border_value=1)
# NOTE: 'border_value' for erosion should be 'True'
# Step 2: Label the eroded pore space
# NOTE: Assume the NW phase reservoir is at the first layer of the domain
# i.e. at seg_image[0,:] - adjust it if this does not suit your need
# For erosion, assume that diagonals are not considered
seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(2,1))
#seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(2,2))
# NOTE: Here I assume the inlet is at the first layer of the array's axis=2 (i.e. domain[0,:,:])\
# You can always change to any other layers as the inlet for this drainage.
label_check = seg_im_ero_label_temp[0,seg_im_ero_label_temp[0,:]!=0]
label_check = np.unique(label_check)
# NOTE the following lines are only for your to check things
# ******************** For check *******************************#
# It assign the labelled array as: NW -> 1, W -> 2, Solid -> 0
#seg_im_ero_label_show = seg_im_ero_label.copy()
#seg_im_ero_label_show[seg_im_ero_label_show !=1] = 2
#seg_im_ero_label_show[np.logical_not(seg_image_2d)]=0
# ******************** End: for check **************************#
seg_im_ero_label = np.zeros_like(seg_im_ero_label_temp,dtype=bool)
for labels in label_check:
seg_im_ero_label = np.logical_or(seg_im_ero_label,seg_im_ero_label_temp==labels)
#seg_im_ero_label = seg_im_ero_label.astype(np.uint8)
# Step 3: perform dilation on the labelled pore space
seg_im_ero_label_dil = morphology.binary_dilation(seg_im_ero_label,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of 'seg_im_ero_label_dil' is 'bool'
seg_im_ero_label_dil[seg_image_input<=0.0]=False
Sw = 1.0 - seg_im_ero_label_dil.sum()/pore_vol
else:
pore_vol = 1.0*(seg_image_input>0.0).sum()
radius = R_critical
print 'Morphological Drainage: processing critical radius: '+str(radius)+' now......'
# Step 1.1: Create structuring element
domain_size = int(np.rint(radius*2)+2)
grid = np.indices((domain_size,domain_size,domain_size))
mk_circle = (grid[0]-domain_size/2)**2 + (grid[1]-domain_size/2)**2 + (grid[2]-domain_size/2)**2 <= radius**2
circle = np.zeros((domain_size,domain_size,domain_size),dtype=np.uint8)
circle[mk_circle]=1
circle = extract_shape(circle).astype(bool)
# Step 1.2: Perform erosion on the pore space
# NOTE: the dtype of 'seg_im_ero' is 'bool'
seg_im_ero = morphology.binary_erosion(seg_image_input>0.0,structure=circle,border_value=1)
# NOTE: 'border_value' for erosion should be 'True'
# Step 2: Label the eroded pore space
# NOTE: Assume the NW phase reservoir is at the first layer of the domain
# i.e. at seg_image[0,:] - adjust it if this does not suit your need
# For erosion, assume that diagonals are not considered
seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(3,1))
#seg_im_ero_label_temp,num_features = measurements.label(seg_im_ero,structure=morphology.generate_binary_structure(3,3))
# NOTE: Here I assume the inlet is at the first layer of the array's axis=2 (i.e. domain[0,:,:])\
# You can always change to any other layers as the inlet for this drainage.
label_check = seg_im_ero_label_temp[0,seg_im_ero_label_temp[0,:]!=0]
label_check = np.unique(label_check)
# NOTE the following lines are only for your to check things
# ******************** For check *******************************#
# It assign the labelled array as: NW -> 1, W -> 2, Solid -> 0
#seg_im_ero_label_show = seg_im_ero_label.copy()
#seg_im_ero_label_show[seg_im_ero_label_show !=1] = 2
#seg_im_ero_label_show[np.logical_not(seg_image_2d)]=0
# ******************** End: for check **************************#
seg_im_ero_label = np.zeros_like(seg_im_ero_label_temp,dtype=bool)
for labels in label_check:
seg_im_ero_label = np.logical_or(seg_im_ero_label,seg_im_ero_label_temp==labels)
#seg_im_ero_label = seg_im_ero_label.astype(np.uint8)
# Step 3: perform dilation on the labelled pore space
seg_im_ero_label_dil = morphology.binary_dilation(seg_im_ero_label,structure=circle,border_value=0)
# NOTE: 'border_value' for dilation should be 'False'
# NOTE: the dtype of 'seg_im_ero_label_dil' is 'bool'
seg_im_ero_label_dil[seg_image_input<=0.0]=False
Sw = 1.0 - seg_im_ero_label_dil.sum()/pore_vol
#end if
return Sw
#end def
def extract_shape(domain):
if domain.ndim == 3:
where_tube = np.where(domain)
zStart = where_tube[0].min()
zEnd = where_tube[0].max()
yStart = where_tube[1].min()
yEnd = where_tube[1].max()
xStart = where_tube[2].min()
xEnd = where_tube[2].max()
domain_seg = domain[zStart:zEnd+1,yStart:yEnd+1,xStart:xEnd+1].copy()
# IMPORTANT: if you have "domain_seg = domain[yStart:yEnd+1,xStart:xEnd+1]"
# then "domain_seg" is only a view of domain, and later on you have
# any changes on your "domain_seg", the "domain" will also be changed
# correspondingly, which might introduce unexpected conflicts and errors
else: # domain.ndim == 2
where_tube = np.where(domain)
yStart = where_tube[0].min()
yEnd = where_tube[0].max()
xStart = where_tube[1].min()
xEnd = where_tube[1].max()
domain_seg = domain[yStart:yEnd+1,xStart:xEnd+1].copy()
return domain_seg
#end def