From 0a769cd6d74a88f27b4afd2bf6a8ecee5a49bc2b Mon Sep 17 00:00:00 2001
From: German <german.martinez-carvajal@irstea.fr>
Date: Tue, 19 Jun 2018 09:52:16 +0200
Subject: [PATCH] 18_06_18
0) Analysis_EM_&_Deposit_Height.py
A new program to analyze the FDEM mapping (in french EM)
FDMI survey: a correlation between deposit height and electric
conductivity mapping and some graphs
It is not possible to observe the correlation, although the relationship
is very intuitive.
1)In PoreSizeDistriution_(ErosionDiltation)_Jython.py
change separator fro, "," into ";" because it is the best way to work in
windows
2)Segmentation_Hashemi.py
I wrote better some explanations
3)Test_Speed_ForLoop_vs_Map.py
Testing speed, trying to understand the performance for vectorized
methods
4)In VER-Aleatoryboxes.py
I changed the source library of the randol function
I added some lines to compute the duration (in seconds) of some parts of
the program
5) The new program VER.py
is intended to be a new version of VER program with a better code, more
comprehensible
---
Analysis_EM_&_Deposit_Height.py | 142 ++++++++++++++++++
...zeDistriution_(ErosionDiltation)_Jython.py | 6 +-
Segmentation_Hashemi.py | 41 ++---
Test_Speed_ForLoop_vs_Map.py | 89 +++++++++++
VER-Aleatoryboxes.py | 67 +++++----
VER.py | 21 +++
6 files changed, 314 insertions(+), 52 deletions(-)
create mode 100644 Analysis_EM_&_Deposit_Height.py
create mode 100644 Test_Speed_ForLoop_vs_Map.py
create mode 100644 VER.py
diff --git a/Analysis_EM_&_Deposit_Height.py b/Analysis_EM_&_Deposit_Height.py
new file mode 100644
index 0000000..fca9f57
--- /dev/null
+++ b/Analysis_EM_&_Deposit_Height.py
@@ -0,0 +1,142 @@
+########################################################
+# Analysis of the relationship between height of deposit
+# layer and electric conductivity (EM)
+########################################################
+
+def Fusion_Frames (Frame_EM, Frame_Deposit):
+ import pandas as pd
+ import numpy as np
+ """ Returns a new frame where a column with the data of "Deposit Heights" belonging to Frame_Deposit
+ is added to Frame_EM at the corresponding coordinates (# of filter, X and Y).
+ If data is missing, the value is filled with np.nan
+
+ Output: Data_EM_DH (pandaDataFrame) "
+ """
+ # FUSIONING BOTH FRAMES
+ # The number of points in both frames is different.
+ # Here I fuse the frames if the coordinates are the same
+ # if not a np.nan is stored
+
+ # Selecting the columns corresponding to the coordinates
+ CoordinatesEM = Frame_EM [["Filtre_#","x_[m]_by_hand","y_[m]_by_hand"]]
+ CoordinatesDP = Frame_Deposit[["Filtre_#","x_[m]_by_hand","y_[m]_by_hand"]]
+ # Initializing the new columns Heights and Comments
+ Dep_Heights = np.zeros(len(Frame_EM)) # initialize Deposit Heights
+ Dep_Heights.fill(np.nan)
+ # In this loop I import the values of Height and comment if coordinates are the same
+ # I need two for loops to study all possibilies
+ for index1 in CoordinatesEM.index:
+ for index2 in CoordinatesDP.index:
+ Same_Coords = CoordinatesEM.ix[index1] == CoordinatesDP.ix[index2] # this is a series-type object
+ Same_Coords = Same_Coords.all() #= True if the three (#Filter, X and Y) coordinates are equal
+ if Same_Coords:
+ Dep_Heights[index1] = Frame_Deposit.ix[index2,["h(cm)" ]]
+ CoordinatesDP = CoordinatesDP.drop(index2)
+ break # break the loop
+ # Checking if all deposit depths were imported
+ if len(CoordinatesDP) != 0:
+ print('WARNING, fuision of data is not complete !')
+ print('There are Deposit Heights data whose coordinates are not present in the EM_Data')
+ else:
+ print('Fusion of data OK!')
+ Data_EM_DH = Frame_EM.copy()
+ Data_EM_DH["h(cm)" ] = Dep_Heights # Deposit Heights
+ return(Data_EM_DH)
+
+def Figure_1 (Frame):
+ import matplotlib.pyplot as plt
+ """ Plots the a scatter of electric conductivity and Deposit Height
+ taking all the points into account """
+ Figure, Ax = plt.subplots()
+ Ax.scatter (Frame["h(cm)" ], Frame["Cond_[mS/m]"])
+ Ax.set_xlabel('Deposit Height (cm)')
+ Ax.set_ylabel('Conductivity (mS/m)')
+ Ax.set_title ('Scatter with all points')
+ plt.show()
+
+def Figure_2 (Frame):
+ import matplotlib.pyplot as plt
+ """ Plots the a scatter of electric conductivity and Deposit Height
+ taking points that were not affected by the presence of a metalic
+ object while performing the EM """
+ # Deleting non consisting values
+ Outliers1 = Frame["Inphase_[ppt]"] > 3 #
+ Outliers2 = Frame["x_[m]_by_hand"] == 0 # Parce qu'il sont proche d'un talud (surface non plate)
+ Outliers3a= Frame["y_[m]_by_hand"] == 0
+ Outliers3b= Frame["Filtre_#"] == 1
+ Outliers3 = Outliers3a & Outliers3b
+ Outliers = Outliers1 | Outliers2 | Outliers3 # Parce qu'il sont proche d'un talud (surface non plate)
+ Frame_Good= Frame[~Outliers]
+ Figure, Ax = plt.subplots()
+ Ax.scatter (Frame_Good["h(cm)" ], Frame_Good["Cond_[mS/m]"])
+ Ax.set_xlabel('Deposit Height (cm)')
+ Ax.set_ylabel('Conductivity (mS/m)')
+ Ax.set_title ('Scatter without interference points')
+ plt.show()
+ return(Frame_Good)# the Frame without the outliers
+
+def Figure_3 (Frame, Filter_num):
+ import matplotlib.pyplot as plt
+ from matplotlib import cm
+ from mpl_toolkits.mplot3d import Axes3D
+ Selection = (Frame["Filtre_#"] == Filter_num)
+ Frame = Frame[Selection]
+ Fig = plt.figure()
+ Ax = Axes3D(Fig)
+ X = Frame["x_[m]_by_hand"]
+ Y = Frame["y_[m]_by_hand"]
+ Z = Frame["Cond_[mS/m]" ]
+ X, Y, Z = create_2D_arrays(X,Y,Z)
+ Ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.viridis)
+ pass
+def create_2D_arrays(X,Y,Z):
+ """Transforms lists of values X,Y Z, into 2D arrays, using interpolation
+ to use them when you want to draw 3D plots or contours,
+ See: https://stackoverflow.com/questions/9008370/python-2d-contour-plot-from-3-lists-x-y-and-rho """
+ import scipy.interpolate
+ import numpy as np
+ # Set up a regular grid of interpolation points
+ xi, yi = np.linspace(X.min(), X.max(), 100), np.linspace(Y.min(), Y.max(), 100)
+ xi, yi = np.meshgrid(xi, yi)
+ # Interpolate
+ rbf = scipy.interpolate.Rbf(X, Y, Z, function='linear')
+ zi = rbf(xi, yi)
+ return (xi, yi, zi)
+################################################################################
+#
+################################################################################
+# Reading File frame with EM data (Carte à Electromagnetisme Multifrequentiel d'un FPR Filtre Planté des Roseaux)
+import pandas as pd
+File_Path = '/home/german.martinez-carvajal/Desktop/2018/EM_Electromag_multifrequentiel/180329_Montromant_EM/donees_EM_correction_tacheo'
+File_Name = File_Path + "/Data_Frame_EM_dots.odf.csv"
+Frame_EM = pd.read_csv(File_Name, sep = ',')
+# Reading Data frame with Deposit Height Data (Data noted by hand in the same FPR )
+File_Path = '/home/german.martinez-carvajal/Desktop/2018/EM_Electromag_multifrequentiel/180329_Montromant_EM/Hauteurs_depot'
+File_Name = File_Path + "/Hauteur_depot.csv"
+Frame_Deposit = pd.read_csv(File_Name, sep = ',')
+
+# Creating the new frame and saving it
+print("Fusing EM and Depotit data")
+Frame_EM_DH = Fusion_Frames (Frame_EM, Frame_Deposit)
+from os import chdir
+File_Path = '/home/german.martinez-carvajal/Desktop/2018/EM_Electromag_multifrequentiel/180329_Montromant_EM/Data_together'
+chdir (File_Path)
+Frame_EM_DH.to_csv('Frame_EM_DH.csv')
+
+# Plotting
+import matplotlib.pyplot as plt
+plt.close("all")
+# Figure 1 - scatter with all points
+Figure_1(Frame_EM_DH)
+# Figure 2 - scatter without interferences
+Frame_EM_DH_Good = Figure_2(Frame_EM_DH)
+# Figure 3 - scatter without interferences
+Figure_3(Frame_EM_DH_Good, Filter_num = 1)
+
+
+
+
+
+
+
+
diff --git a/PoreSizeDistriution_(ErosionDiltation)_Jython.py b/PoreSizeDistriution_(ErosionDiltation)_Jython.py
index eab1ea3..af9087c 100644
--- a/PoreSizeDistriution_(ErosionDiltation)_Jython.py
+++ b/PoreSizeDistriution_(ErosionDiltation)_Jython.py
@@ -66,9 +66,9 @@ def Pore_Size_Distribution(Dir, FileName):
# Saving
chdir(Dir)
File = open('Cumul_distribution.csv', mode = 'w')
- File.write('Radius (vox)' + ',' + 'Cumul Volume (vox)' + '\n')
+ File.write('Radius (vox)' + ';' + 'Cumul Volume (vox)' + '\n')
for i in range(len(Radii)):
- File.write(str(Radii[i]) + ',' + str(Cum_Volumes[i]) + '\n')
+ File.write(str(Radii[i]) + ';' + str(Cum_Volumes[i]) + '\n')
File.close()
print("Finished !")
@@ -80,9 +80,7 @@ import datetime
Start = datetime.datetime.today()
# Directory and File Name
Dir = '/home/german.martinez-carvajal/Desktop/2018/Pore_Size_Distribution/Test'
-#Dir = '/home/german.martinez-carvajal/Desktop/2017/Tomographie/170411COP/Segmentation(Hashemi)'
FileName = Dir + "/" + "Binary.tif"
-#FileName = Dir + "/" + "Air.tif"
# Calling the function
Pore_Size_Distribution(Dir, FileName)
# Report end-time of the execution
diff --git a/Segmentation_Hashemi.py b/Segmentation_Hashemi.py
index a6980e6..cea2d2f 100644
--- a/Segmentation_Hashemi.py
+++ b/Segmentation_Hashemi.py
@@ -1,26 +1,29 @@
from __future__ import division
-def Segmentation_Hashemi(Dir, Scan_File,PVE_File, Mu = None):
- """ This function creates a 3-phases, segmentation from a 3D computed tomography
- Those 3 phases are often: Air, Organic Matter and Gravel
- Based on: A tomographic imagery segmentation methodology for three-phase geomaterials based on simultaneous region growing
- October 2013Acta Geotechnica 9(5)
-
- I takes 4 inputs
- Dir : The directory were the files are read and written
- Scan_File : .tif 3D-image in 8-bit or 16-bit. It should be filtered before (Use 3D median filter in ImageJ)
- PVE_File : .tif 3D-image in 8-bit (same size as Scan_File). It is binary, white pixels (255 in 8-bit) represent pixels of the Partial Volume Effect
- In imageJ Use the command: Process-Find Edges (Warning, this is a 2D operation)
+def Segmentation_Hashemi(Dir, Scan_File, PVE_File, Mu = None):
+ """ This function creates a 3-phase-segmentation from an image issue of 3D computed tomography
+ Those 3 phases are often: Air, Organic Matter (OM) and Gravel
+
+ The algorithm is based on: A tomographic imagery segmentation methodology for three-phase geomaterials based on simultaneous region growing
+ October 2013Acta Geotechnica 9(5) Hashemi et al
+
+ The funcion takes 4 inputs
+ Dir : The directory where the files Scan_File and PVE_File are read and where results are written
+ Scan_File : .tif 3D-image in 8-bit or 16-bit. It is a previously filtered 3D Image (Use 3D median filter in ImageJ)
+ PVE_File : .tif 3D-image in 8-bit (same size as Scan_File). It is a binary 3D image,
+ white pixels (i.e. value 255 in 8-bit) represent pixels of the Partial Volume Effect
+ To obtain the a PVE Image: In ImageJ use the command: Process-Find Edges (Warning, this is a 2D operation)
Then binarize using the defalut ImageJ thresholding method
- Mu : a list with three integers values that will be the initialisation in the fitting of gaussians of the gray value histogram , if None, they
- will be initialised by some values depending our experience
+ Mu : a list with three integers values that will be the initialisation
+ in the fitting of gaussians of the gray value histogram of ScanFile, if Mu == None, they
+ will be initialised by some values depending on our experience"
Outputs
- 5 figures with histograms to understand how the seed image was built and a track of the growing phase
- The seed Image (following Hashemi's method)
- An image for every 5th growing steps
- A inmage for interface pixels attribution
- The final segmentation (labels are 0 for air, 128 for OM and 255 for Gravel)
- Some messages about the performance of the segmentation
+ 5 figures with histograms to understand how the seed-image was built and a track of the growing phase
+ -The seed Image in .tif format (following Hashemi's method)
+ -An image for every 5th growing steps in .tif format
+ -An image for interface pixels attribution in .tif frmat
+ -THE MOST IMPORTANT: The final segmentation (WHOSE labels are 0 for air, 128 for OM and 255 for Gravel) in .tif format
+ -Some messages about the performance of the segmentation in a .csv file
"""
from os import chdir, path
diff --git a/Test_Speed_ForLoop_vs_Map.py b/Test_Speed_ForLoop_vs_Map.py
new file mode 100644
index 0000000..0f68619
--- /dev/null
+++ b/Test_Speed_ForLoop_vs_Map.py
@@ -0,0 +1,89 @@
+import numpy as np
+import pandas as pd
+from numpy.random import randint as randint
+import datetime
+
+number_radii = 20
+number_boxes = 20
+size = (10,10)
+
+def create_list_of_boxes():
+ print('Creating list of voxels')
+ global list_of_boxes
+ start = datetime.datetime.today()
+ A = randint(0,3, size = size, dtype = np.uint8)
+ list_of_boxes = [[A for j in range (number_boxes)] for i in range (number_radii)]
+ end = datetime.datetime.today()
+ duration = end - start
+ duration = str(duration)
+ print('Boxes CreationV1 - as a list of lists')
+ print('(microseconds): ', duration)
+
+def create_list_of_boxesV2():
+ global Frame
+ start = datetime.datetime.today()
+ Radii = []
+ Box_Number = []
+ Boxes = []
+ A = randint(0,3, size = size, dtype = np.uint8)
+ for Radius in range(number_radii):
+ for Box in range(number_boxes):
+ Radii .append(Radius)
+ Box_Number.append(Box)
+ Boxes .append(A)
+ Frame = pd.DataFrame({'Radii':Radii, 'Box_Number':Box_Number,'Boxes':Boxes})
+ end = datetime.datetime.today()
+ duration = end - start
+ duration = str(duration)
+ print('Boxes CreationV2 - as a panda frame passing through a for')
+ print('(microseconds): ', duration)
+
+
+def Compute_Volume_Fractions (box):
+ assert (type(Box) == np.ndarray), 'Input "Box" must be a numpy.ndarray'
+ Air_Voxels = np.count_nonzero((box==0)) #Air
+ OrM_Voxels = np.count_nonzero((box==1)) #OrM
+ Grv_Voxels = np.count_nonzero((box==2)) #Grv
+ Test = ( np.product(box.shape) == ( Air_Voxels + OrM_Voxels + Grv_Voxels ) )
+ if not (Test == True):
+ raise AssertionError("You must verify that all phase-labels in Box are 0, or, 128, or 255")
+ Volume = Air_Voxels + OrM_Voxels + Grv_Voxels
+ VF_Air, VF_OrM, VF_Grv = Air_Voxels/Volume, OrM_Voxels/Volume, Grv_Voxels/Volume
+ return (VF_Air, VF_OrM, VF_Grv)
+
+# Creating list of boxes
+# First way
+create_list_of_boxes()
+# Second way
+create_list_of_boxesV2()
+
+print('Computing Volume Fractions: ')
+# Way 1
+print ('Way 1: For_loop (microseconds): ' , str(duration))
+start = datetime.datetime.today()
+VF = [[Compute_Volume_Fractions(Box) for Box in list_of_boxes [i]] for i in range (number_radii)]
+VF = np.array(VF)
+end = datetime.datetime.today()
+duration = end - start
+duration = str(duration)
+
+# Way 2
+print ('Way 2: For + Map (microseconds): ' , str(duration))
+start = datetime.datetime.today()
+VF2 = [list(map(Compute_Volume_Fractions, list_of_boxes[i])) for i in range (number_radii) ]
+VF2 = np.array(VF2)
+end = datetime.datetime.today()
+duration = end - start
+duration = str(duration)
+
+# Way 3
+print ('Way 3: Map Only (microseconds): ' , str(duration))
+start = datetime.datetime.today()
+VF3 = list(map(Compute_Volume_Fractions, Frame['Boxes']))
+VF3 = np.array(VF3)
+end = datetime.datetime.today()
+duration = end - start
+duration = str(duration)
+
+
+print ('ok?: ', np.alltrue(VF ==VF2))
\ No newline at end of file
diff --git a/VER-Aleatoryboxes.py b/VER-Aleatoryboxes.py
index 2c74172..e0c869a 100644
--- a/VER-Aleatoryboxes.py
+++ b/VER-Aleatoryboxes.py
@@ -3,23 +3,17 @@ from os import chdir
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
-from numpy.random import random_integers as rd_int
-plt.close('all')
+from numpy.random import randint as rd_int
+import datetime
-Directory_reading = "/home/german.martinez-carvajal/Documents/Tomographie/170418ALL/ALL02/Small Sample/Growing-Hashemi"
-Directory_saving = "/home/german.martinez-carvajal/Documents/Tomographie/170418ALL/ALL02/Small Sample/Growing-Hashemi"
-Image_Name = 'Hashemi_Segmentation'+".tif"
-chdir(Directory_reading)
-Segmentation = tifffile.imread (Image_Name)
def Compute_Volume_Fractions (box):
Air_Voxels = np.count_nonzero((box==0 )) #Air
- OrM_Voxels = np.count_nonzero((box==128)) #OrM
+ OrM_Voxels = np.count_nonzero((box==128)) #OrM or OM (Organic Matter)
Grv_Voxels = np.count_nonzero((box==255)) #Grv
Test = ( np.product(box.shape) == ( Air_Voxels + OrM_Voxels + Grv_Voxels ) )
if not (Test == True):
raise AssertionError("You must verify that all phase-labels in Box are 0, or, 128, or 255")
-
Volume = Air_Voxels + OrM_Voxels + Grv_Voxels
VF_Air, VF_OrM, VF_Grv = Air_Voxels/Volume, OrM_Voxels/Volume, Grv_Voxels/Volume
return (VF_Air, VF_OrM, VF_Grv)
@@ -77,7 +71,6 @@ def Create_2D_VER_Graph_AB (Segmentation, dsave, Z_slice, Resolution, N_box = 10
Center = (int(0.5*Segmentation.shape[0]),int(Segmentation.shape[1]*0.5))
# List of radii of the aleatory boxes
- print("Creating Radii")
global Radii
Points_number = Radii_number # Points number in the graph (same as number of radii)
Radiux_Max = int(0.5*np.min(Segmentation.shape))
@@ -146,16 +139,31 @@ def Create_2D_VER_Graph_AB (Segmentation, dsave, Z_slice, Resolution, N_box = 10
Bigger = Make_Simetric_Image_2D (Segmentation)
# Sampling the of boxes
Boxes2D = [[Bigger[X_left[i][j]:X_right[i][j], Y_left[i][j]:Y_right[i][j]] for j in range(Number_of_Boxes) ] for i in range (Radii.size)]
-# # Adding a las element to be sure that the last group of boxex is the whole 2D image
-# Radii = np.append(Radii,0.5*np.min(Segmentation.shape))
-# Boxes2D.append([Segmentation]*Number_of_Boxes)
-
- # Computing volume fractions
- print("Computing volume fractions")
+ ## Adding a las element to be sure that the last group of boxex is the whole 2D image
+ #Radii = np.append(Radii,0.5*np.min(Segmentation.shape))
+ #Boxes2D.append([Segmentation]*Number_of_Boxes)
+
+ #Boxes2D = np.array(Boxes2D)
+ #Computing volume fractions
+ print("Computing volume fractions")
+ start = datetime.datetime.today()
global Means_Air, Means, Std, P25, P50, P75, Box_lenghts_Rel, VF
VF = [[Compute_Volume_Fractions(Box) for Box in Boxes2D [i]] for i in range (Radii.size)]
VF = np.array(VF)
+ end = datetime.datetime.today()
+ duration = end - start
+ duration = str(duration)
+ print ('Duration for double loop: ' , str(duration))
+
+ start = datetime.datetime.today()
+ VF2 = [list(map(Compute_Volume_Fractions, Boxes2D [i])) for i in range (Radii.size) ]
+ VF2 = np.array(VF2)
+ end = datetime.datetime.today()
+ duration = end - start
+ duration = str(duration)
+ print ('Duration for loop + map : ' , str(duration))
+ print ('ok?: ', np.alltrue(VF ==VF2))
# Computing means and standard deviations
print("Computing statistics")
@@ -394,8 +402,8 @@ def Repeat_for_many_slices(Start_Slice, Segmentation, Num_slices, Dir_saving, Re
# Closing preexisting figures
plt.close('all')
# Directories
-Dir_reading = "/home/german.martinez-carvajal/Documents/Tomographie/170411COP/Segmentation(Hashemi)"
-Dir_saving = "/home/german.martinez-carvajal/Documents/Tomographie/170411COP/VER"
+Dir_reading = "/home/german.martinez-carvajal/Desktop/2017/Tomographie/170411COP/Segmentation(Hashemi)"
+Dir_saving = "/home/german.martinez-carvajal/Desktop/2017/Tomographie/170411COP/VER"
# Reading segmented image
chdir(Dir_reading)
Image_Name = 'Sample-100slices-Values(0-128-255)'+".tif"
@@ -407,17 +415,17 @@ Z_slice = 0
## Call function for the test of the number of aleatory boxes
#Means_Air, Errors = Study_of_number_of_AB (Segmentation, Dir_saving, Resolution, 'Reflect')
-## Call function AB : Aleatory Boxes
-#Create_2D_VER_Graph_AB (Segmentation[Z_slice], Dir_saving, Z_slice, Resolution, N_box = 2000, check_centers = False, Save_Graphs = True, Radii_number = 20, Mode = "Periodic")
+# Call function AB : Aleatory Boxes
+Create_2D_VER_Graph_AB (Segmentation[Z_slice], Dir_saving, Z_slice, Resolution, N_box = 200, check_centers = False, Save_Graphs = True, Radii_number = 20, Mode = "Periodic")
-# Call function to repeat VER evaluation for several slices
-chdir(Dir_reading)
-Image_Name = 'HashemiConcatenated(0-128-255)'+".tif"
-Segmentation = tifffile.imread (Image_Name)
-Num_Slices = 4
-Dir_saving = "/home/german.martinez-carvajal/Documents/Tomographie/170411COP/VER/Depth_Analysis"
-Start_Slice = 400
-Repeat_for_many_slices(Start_Slice, Segmentation, Num_Slices, Dir_saving, Resolution, N_box = 2000, Mode = 'Reflect' )
+## Call function to repeat VER evaluation for several slices
+#chdir(Dir_reading)
+#Image_Name = 'HashemiConcatenated(0-128-255)'+".tif"
+#Segmentation = tifffile.imread (Image_Name)
+#Num_Slices = 4
+#Dir_saving = "/home/german.martinez-carvajal/Documents/Tomographie/170411COP/VER/Depth_Analysis"
+#Start_Slice = 400
+#Repeat_for_many_slices(Start_Slice, Segmentation, Num_Slices, Dir_saving, Resolution, N_box = 2000, Mode = 'Reflect' )
## Call function for Montromant station
## Directories
@@ -433,4 +441,5 @@ Repeat_for_many_slices(Start_Slice, Segmentation, Num_Slices, Dir_saving, Resolu
# Create_2D_VER_Graph_AB (Segmentation[Z_slice], Dir_saving, Z_slice, Resolution, N_box = 2000, check_centers = False, Save_Graphs = True, Radii_number = 20, Mode = Mode)
#plt.close('all')
-print("FINISHED")
\ No newline at end of file
+print("FINISHED")
+plt.close("all")
\ No newline at end of file
diff --git a/VER.py b/VER.py
new file mode 100644
index 0000000..243a270
--- /dev/null
+++ b/VER.py
@@ -0,0 +1,21 @@
+from __future__ import division
+
+def Compute_Volume_Fractions (box):
+ """ Returns a tupple with 3 the value of three volume fractions (VF_Air, VF_OrM, VF_Grv)
+ from a input "box" (type = numpy.ndarray, dtype = numpy.uint8) whose labels are EXCLUSIVELY
+ 0 for phase1 (Air), 128 for phase2 (OrM organic matter), 255 for phase3 (Gravel) """
+ import numpy as np
+ # Checking type of "box"
+ if not (type(box) == np.ndarray):
+ raise TypeError("Input MUST be a numpy.ndarray")
+ # Counting pixels belonging to each phase
+ Air_Voxels = np.count_nonzero((box==0 )) #Air
+ OrM_Voxels = np.count_nonzero((box==128)) #OrM
+ Grv_Voxels = np.count_nonzero((box==255)) #Grv
+ Test = (np.product(box.shape) == (Air_Voxels + OrM_Voxels + Grv_Voxels))
+ if not (Test == True):
+ raise AssertionError("You must verify that all phase-labels in Box are 0, or, 128, or 255")
+ Volume = Air_Voxels + OrM_Voxels + Grv_Voxels
+ VF_Air, VF_OrM, VF_Grv = Air_Voxels/Volume, OrM_Voxels/Volume, Grv_Voxels/Volume
+ return (VF_Air, VF_OrM, VF_Grv)
+
--
GitLab