diff --git a/Analysis_EM_&_Deposit_Height.py b/Analysis_EM_&_Deposit_Height.py new file mode 100644 index 0000000000000000000000000000000000000000..fca9f576d4f2676e7b7062283b29eec3bd1b0f0e --- /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 eab1ea3953d693d40fe5c364f01921eb0786bdc5..af9087c442fb133e373936de9f465f63ccf83709 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 a6980e609f1fa668674ce6d6d0005c0f8d6e3217..cea2d2fe23a7274d4f477e151743f13ecc03465b 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 0000000000000000000000000000000000000000..0f686199929360188a4efc02eed725377fe55d8c --- /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 2c741725158f302cf258c66e59cc5102c099ee74..e0c869ad43aed1b351b75b0d641045aec870ee14 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 0000000000000000000000000000000000000000..243a27091362cc75f90b2888bac7d8ea409648bf --- /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) +