Commit 6c6767a3 authored by Martinez-Carvajal German's avatar Martinez-Carvajal German
Browse files

adding all possible combinations to compute

specific surfaces where the FM/voids surface
is divided by
1) FM volume
2) Air Volume
3) FM+Air Volume
4) Air + Grv + FM volume
parent 91d0bb70
......@@ -41,7 +41,9 @@ def SurfaceCouchyCrofton_2(Image, Resolution, p2, Surface_of_p1):
global Progress
assert(Image.shape == Surface_of_p1.shape)
assert(p2 != 0), 'arbitrarily, the value of phase 2 (p2) MUST be different from 0, CHANGE the labels from your original image'
# Relabeling
num_sur_voxels = np.count_nonzero(Surface_of_p1 == 255)
......@@ -240,9 +242,11 @@ def plot_ssp(path, file):
data = pd.read_csv(file, sep = ';' )
plt.close('all')
figure, ax = plt.subplots()
ax.plot(data['SS (mm-1) FM'], data['Depth(mm)'], linewidth = 1, linestyle = '-', color = 'k', label = "per volume of FM")
ax.plot(data['SS (mm-1) Bulk'], data['Depth(mm)'], linewidth = 1, linestyle = '--', color = 'k', label = "per bulk volume")
ax.plot(data['SS (mm-1) / FM'], data['Depth(mm)'], linewidth = 1, linestyle = '-', color = 'red', label = "per FM volume")
ax.plot(data['SS (mm-1) / Bulk'], data['Depth(mm)'], linewidth = 1, linestyle = '-', color = 'purple', label = "per bulk volume")
ax.plot(data['SS (mm-1) / Voids'], data['Depth(mm)'], linewidth = 1, linestyle = '-', color = 'blue', label = "per voids volume")
ax.plot(data['SS (mm-1) / Voids + FM'], data['Depth(mm)'], linewidth = 1, linestyle = '-', color = 'black', label = "per voids + FM volume")
ax.set_xlim(0)
ax.legend(loc = 'best')
plt.xlabel(r'Specific Surface $mm^{2}/mm^{3}$')
......@@ -288,7 +292,7 @@ def compute_specific_surface_profile(start_VS, path_input, path_save, file_name,
list_messages_out = add_message('Reading succesful', list_messages_out, path_save, file_name)
# extracting surface pixels from phase1
# extracting surface pixels from phase1 (p1)
if flag_SV == True:
print("Extracting surface voxels...")
SurfaceVoxels = extract_surface_voxels(image, p1, file_name)
......@@ -299,73 +303,105 @@ def compute_specific_surface_profile(start_VS, path_input, path_save, file_name,
# Relabel (to avoid pixels outside the cylinder)
if flag_relabel == True:
image = relabel_the_outside(image, label_outside, rw)
#chdir(path_save)
#tifffile.imsave('relabeled_{}'.format(file_name), image)
# computing indexes for each vertical section
# computing indexes for each vertical section
lenght = len(image)
limits_of_VS = np.arange(0,lenght+size_VS,size_VS) # Limits stores the index of the boxes where a profile will be performed
# computing the profile
SpcSrfcs = list() # Voids-FM surface / FM volume
# Initilializing the variables that will be computed for each section
# Volumes
VolumesOrM = list()
VolumesAir = list()
Volumes = list() # sum of air, FM and gravel
# Surface
SurfAir_OrM= list()
VolumesOrM = list()
# Specific surfaces (depending "volume" you will use to divide)
SpcSrfcs_1 = list() # Voids-FM surface / FM volume
SpcSrfcs_2 = list() # Voids-FM surface / air volume
SpcSrfcs_3 = list() # Voids-FM surface / volume air and FM volume
SpcSrfcs_4 = list() # Voids-FM surface / volume air, FM, and gravel volume
# Depth
Depth = list()
# Count of voxels located at the surface
num_sur_voxs = list()
# Contributions of each family in the couchy crofton's method
F1s = list()
F2s = list()
F3s = list()
Volumes = list()
SpcSrfcs_2 = list() # Voids-FM surface / bulk volume
# Computing the specific surface small slices with size = SizeBox
# Computing the specific surface small sections with thickness = SizeBox
for i in range(start_VS,len(limits_of_VS)-1):
# Saving messages
message = 'Computing profile, step {}/{}'.format(i+1,len(limits_of_VS)-1 )
list_messages_out = add_message(message, list_messages_out, path_save, file_name)
print('Step: {}/{}'.format(i+1,len(limits_of_VS)-1 ))
# Moving the window downwars to perform the computations
Selection1 = image [limits_of_VS[i]:limits_of_VS[i+1]].copy()
Selection2 = SurfaceVoxels[limits_of_VS[i]:limits_of_VS[i+1]].copy()
#tifffile.imsave('S1_{}'.format(i), Selection1)
#tifffile.imsave('S2_{}'.format(i), Selection2)
# Volume of Fouling Material in the section
VolumeOrM = np.count_nonzero(Selection1 == p2)
VolumeOrM = VolumeOrM*resolution**3 # in mm3
print('FM Volume of section =', VolumeOrM)
# Volume of air in the section
VolumeAir = np.count_nonzero(Selection1 == p1)
VolumeAir = VolumeAir*resolution**3 # in mm3
# volume of the section (take out pixels outside the cylinder)
Volume = np.count_nonzero(Selection1 != label_outside ) * resolution**3
Volume1 = (np.count_nonzero(Selection1 == 0) + np.count_nonzero(Selection1 == 128) + np.count_nonzero(Selection1 == 255))*resolution**3
assert Volume == Volume1
print('FM Volume of section =', VolumeOrM)
Surface, num_sur_vox, (F1, F2, F3) = SurfaceCouchyCrofton_2(Selection1, resolution, p2, Selection2)
print('Surface of section =', Surface)
print('Surface of section (mm2)=', Surface)
# Storing data
SpecSur = Surface/VolumeOrM # Voids-FM surface / FM volume
SpcSrfcs .append(SpecSur )
SurfAir_OrM.append(Surface )
VolumesOrM .append(VolumeOrM)
# Volumes
VolumesOrM.append(VolumeOrM)
VolumesAir.append(VolumeAir)
Volumes.append(Volume)
SpcSrfcs_2.append(Surface/Volume) # Voids-FM surface / bulk volume
# Suface
SurfAir_OrM.append(Surface)
# Specific surfaces
SpcSrfcs_1.append(Surface/VolumeOrM)
SpcSrfcs_2.append(Surface/VolumeAir)
SpcSrfcs_3.append(Surface/(VolumeOrM+VolumeAir))
SpcSrfcs_4.append(Surface/(Volume))
# Count of surface voxels and contributions
num_sur_voxs.append(num_sur_vox)
F1s.append(F1)
F1s.append(F1)
F2s.append(F2)
F3s.append(F3)
Depth.append(limits_of_VS[i]*(-1)*(resolution)) # in mm)
# Depth
Depth.append(limits_of_VS[i]*(-1)*(resolution)) # in mm
print('value SpecSurf',limits_of_VS[i],'-',limits_of_VS[i+1], '=', SpecSur)
print('value SpecSurf / FM volume',limits_of_VS[i],'-',limits_of_VS[i+1], '=', Surface/VolumeOrM)
# Save and plot (little by little)
chdir(path_save)
Frame = pd.DataFrame()
Frame['Depth(mm)'] = Depth
Frame['SS (mm-1) FM'] = SpcSrfcs
Frame['SS (mm-1) Bulk'] = SpcSrfcs_2
Frame['Surf Voids/OrM (mm^2)']= SurfAir_OrM
Frame['Volumes Bulk (mm^3)'] = Volume
Frame['Volumes OrM (mm^3)'] = VolumesOrM
Frame['Volumes FM (mm^3)'] = VolumesOrM
Frame['Volumes Voids (mm^3)'] = VolumesAir
print(VolumesAir, VolumesOrM)
Frame['Volumes Voids + FM (mm^3)'] = np.array(VolumesAir) + np.array(VolumesOrM)
Frame['Volumes Bulk (mm^3)'] = Volumes
Frame['SS (mm-1) / FM'] = SpcSrfcs_1
Frame['SS (mm-1) / Voids'] = SpcSrfcs_2
Frame['SS (mm-1) / Voids + FM'] = SpcSrfcs_3
Frame['SS (mm-1) / Bulk'] = SpcSrfcs_4
Frame['Count_Sur_Vox'] = num_sur_voxs
Frame['F1'] = F1s
Frame['F2'] = F2s
Frame['F3'] = F3s
Frame['Sur_Vox'] = num_sur_voxs
Name = 'SS_Profile_of_{}_rw_{}.csv'.format(sample_name, rw)
Frame.to_csv(Name, sep = ';')
......@@ -397,17 +433,17 @@ def main():
# Image resolution
resolution = 0.035 #mm per vox
# Labels of the phases
# Labels of the phases p1 is ALWAYS the label for air, and p2 is ALWAYS the label for OrM (Organig matter or Fouling Material)
label_outside = 1 # 1 = outside the cylinder, 0 = voids, 128 = fouling material, 255 = gravel
p1 = 0 # computation will be done for interface between phases p1 and p2
p2 = 128
# Relabel ? (see relabel function)
flag_relabel = True
rw = 1
rw = 1.0
# size_VS
size_VS = 2 # pixels (the same lenght of the smoothing for volume fraction profiles)
size_VS = 3 # pixels (the same lenght of the smoothing for volume fraction profiles)
# start at vertical section #...
# useful if you want to resume the computation from a certain point
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment