Commit 8d4a90a9 authored by Martinez-Carvajal German's avatar Martinez-Carvajal German
Browse files

one file replaces the other

the new one computes everything, global volume fractions, global volumes, and depths
the 3 of them for the deposit layer and the gravel layer once the limits have been established by the vertical volume fraction profiles
parent d3ff1919
# -*- coding: utf-8 -*-
"""
Created on Wed May 22 13:08:22 2019
@author: german.martinez-carvajal
"""
import tifffile
import numpy as np
from os import chdir
import pandas as pd
def sector_mask(shape,centre,radius,angle_range):
"""
CODE FROM STACK OVER FLOW
Return a boolean mask for a circular sector. The start/stop angles in
`angle_range` should be given in clockwise order.
"""
x,y = np.ogrid[:shape[0],:shape[1]]
cx,cy = centre
tmin,tmax = np.deg2rad(angle_range)
# ensure stop angle > start angle
if tmax < tmin:
tmax += 2*np.pi
# convert cartesian --> polar coordinates
r2 = (x-cx)*(x-cx) + (y-cy)*(y-cy)
theta = np.arctan2(x-cx,y-cy) - tmin
# wrap angles between 0 and 2*pi
theta %= (2*np.pi)
# circular mask
circmask = r2 <= radius*radius
# angular mask
anglemask = theta <= (tmax-tmin)
return circmask*anglemask
def relabel_the_outside(image,labels):
""" if pixels outside the cylinder that contains the sample have the same value for air
they will be relabeled with the value label[1]
so, if you follow the conventions for label values
label [0]= 1 = pixels outside the cylinder,
label [1]= 0 = air (i.e voids) pixels
label [2]= 128 = fouling material
label [3]= 255 = gravel
"""
assert (np.count_nonzero(image == 1) == 0 ), 'You must reserve label = 1, for pixels outside the cylinder'
print("relabeling from 0,128,255 to 0,1,128, 255")
# CREATING CIRCULAR MASK
# The cylinder is supposed to be the biggest "circle" inscribed in a squared image
z,x,y = np.shape(image)[0], np.shape(image)[1], np.shape(image)[2]
center = int(x/2), int(y/2)
radius = int((x/2))
#Selecting circular ROI (region of interest)
mask_2D = sector_mask((x,y),center,radius,(0,360 ))
mask_3D = np.zeros(image.shape, dtype = np.bool)
for j in range(z):
mask_3D[j] = mask_2D
image[np.logical_not(mask_3D)] = labels[0] # outside the cilinder label = 1
print("{} pixels have been attributed to tha phase : outside the cylinder".format(np.count_nonzero(image == 1)))
return image
def compute_global_fractions_volumes_depths(r, image, sample_name, start_dl, end_dl, end_gl, path_save):
""" r = resolution in mm/vox
start_dl = int, starting slice of the deposit layer
end_dl = int, ending slice of the deposit layer, start of the gravel layer
end_gl = int, ending of the gravel layer # should be the ending slice of the image
"""
print('computing...')
zones = ['Total','Deposit Layer', 'Gravel Layer']
phases = ['Outside', 'Voids', 'Fouling Material' , 'Gravel']
labels = [1,0,128,255] # each phase has its label
df = pd.DataFrame()
df["Zone"] = zones
depths= []
depths.append((len(image)-start_dl)*r)
depths.append((end_dl-start_dl)*r)
depths.append((len(image)-end_dl)*r)
df["depth (mm)"]= depths
# total volumes in vox (it is necessary to subtract voxels corresponding to the phase outside the cylinder)
global_volume = np.product(image[start_dl:].shape) - np.count_nonzero(image[start_dl:] == 1)
volume_dep_layer = np.product(image[start_dl:end_dl].shape) - np.count_nonzero(image[start_dl:end_dl] == 1)
volume_grv_layer = np.product(image[end_dl:end_gl].shape)- np.count_nonzero(image[end_dl:end_gl] == 1)
# loop over all the phases ('Outside', 'Voids', 'Fouling Material' , 'Gravel')
for i in range(4):
phase = phases[i]
label = labels[i]
fractions = []
volumes = []
if phase == 'Outside':
# global volume fraction
fractions.append(np.count_nonzero(image[start_dl:] == label)/np.product(image[start_dl:].shape))
# volume fraction in the deposit layer
fractions.append(np.count_nonzero(image[start_dl:end_dl] == label)/np.product(image[start_dl:end_dl].shape))
# volume fraction in the gravel layer
fractions.append(np.count_nonzero(image[end_dl:end_gl] == label)/np.product(image[end_dl:end_gl].shape))
else:
# global volume fraction
fractions.append(np.count_nonzero(image[start_dl:] == label)/global_volume)
# volume fraction in the deposit layer
fractions.append(np.count_nonzero(image[start_dl:end_dl] == label)/volume_dep_layer)
# volume fraction in the gravel layer
fractions.append(np.count_nonzero(image[end_dl:end_gl] == label)/volume_grv_layer)
# global volume # IN cm3
volumes.append(np.count_nonzero(image[start_dl:] == label)*((r**3))/1000)
# volume of the deposit layer
volumes.append(np.count_nonzero(image[start_dl:end_dl] == label)*((r**3))/1000)
# volume of the gravel layer
volumes.append(np.count_nonzero(image[end_dl:end_gl] == label)*((r**3))/1000)
df["volume fraction of "+ phase] = fractions
df["volume (cm3) of " + phase ] = volumes
chdir(path_save)
print("Fractions for sample {}\n".format(sample_name), df)
df.to_csv('Global_Volume_Fractions_of_sample_{}.csv'.format(sample_name), sep = ';' )
return (df)
def main():
path_read = '/home/german.martinez-carvajal/Bureau/These/Hashemi_segmentation/Tests/Test9_Rapide/Results/heavy_results'
chdir(path_read) # just to check it exists
path_save = '/home/german.martinez-carvajal/Bureau/These/Total_volumes_&_Global_fractions/test'
chdir(path_save) # just to chech it exists
file_name = 'Hash_Segm_of_Image_Filtered_PVE_50_Tol_0.75_NumVS_1.tif'
sample_name = "Test"
slice_start_deposit_layer = 1
slice_end_deposit_layer = 5
r = 0.035 # mm/vox
# Computing global fractions
chdir(path_read)
print('Reading...', path_read, file_name)
image = tifffile.imread(file_name)
print('Relabeling the outside')
labels = [1,0,128,255] # 1 = outside the cylinder, 0 = voids, 128 = fouling material, 255 = gravel
image = relabel_the_outside(image, labels)
compute_global_fractions_volumes_depths(r, image, sample_name, slice_start_deposit_layer, slice_end_deposit_layer, len(image), path_save)
main()
print('ok global fractions')
# -*- coding: utf-8 -*-
"""
Created on Wed May 22 13:08:22 2019
@author: german.martinez-carvajal
"""
import tifffile
import numpy as np
from os import chdir
import pandas as pd
def compute_volumes(image, labels, resolution):
# resolution in mm/vox
volume_voxels = np.zeros(len(labels))
volume_cm3 = np.zeros(len(labels))
for i in range(len(labels)):
volume_voxels[i] = np.count_nonzero((image==labels[i])) # outside
test = np.product(image.shape) == np.sum(volume_voxels)
if test == False:
raise AssertionError("You must verify that there are EXACTLY {} labels for phases in your image".format(len(labels)))
volume_cm3 = volume_voxels*(resolution**3)/1000
return(volume_cm3)
def compute_global_fractions(image, sample_name, start_dl, end_dl, end_gl, path_save):
""" start_dl = starting slice of the deposit layer
end_dl = ending slice of the deposit layer, start of the gravel layer
end_gl = ending of the gravel layer
"""
print('computing...')
zones = ['Total','Deposit Layer', 'Gravel Layer']
phases = ['Outside', 'Voids', 'Fouling Material' , 'Gravel']
labels = [1,0,128,255] # each phase has its label
df = pd.DataFrame()
df["Zone"] = zones
# total volumes in vox (it is necessary to subtract voxels corresponding to the phase outside the cylinder)
global_volume = np.product(image.shape) - np.count_nonzero(image == 1)
volume_dep_layer = np.product(image[start_dl:end_dl].shape) - np.count_nonzero(image[start_dl:end_dl] == 1)
volume_grv_layer = np.product(image[end_dl:end_gl].shape)- np.count_nonzero(image[end_dl:end_gl] == 1)
# loop over all the phases ('Outside', 'Voids', 'Fouling Material' , 'Gravel')
for i in range(4):
phase = phases[i]
label = labels[i]
fractions = []
# global volume fraction
fractions.append(np.count_nonzero(image == label)/global_volume)
# volume fraction in the deposit
fractions.append(np.count_nonzero(image[start_dl:end_dl] == label)/volume_dep_layer)
# volume fraction in the deposit
fractions.append(np.count_nonzero(image[end_dl:end_gl] == label)/volume_grv_layer)
df[phase] = fractions
chdir(path_save)
print("Fractions for sample {}\n".format(sample_name), df)
df.to_csv('Global_Volume_Fractions_of_sample_{}.csv'.format(sample_name), sep = ';' )
return (df)
def main():
samples = ['MON1A', 'MON1B','MON2A','MON2B','MON3A','MON3B']
paths = {}
for sample in samples:
paths[sample] = '/home/german.martinez-carvajal/Desktop/These/Hashemi_segmentation/Tomos_2018/Manip_juillet/{}_20180723/heavy_results'.format(sample)
files = {}
files['MON1A'] = 'Hashemi_Segmentation_of_MON1A_20180723_ROI_Filtered_aniso_Num_VS_1relabeled.tif'
files['MON1B'] = 'Hashemi_Segmentation_of_MON1B_20180723_ROI_Filtered_aniso_Num_VS_1relabeled.tif'
files['MON2A'] = 'Hashemi_Segmentation_of_MON2A_20180723_ROI_Filtered_aniso_Num_VS_1relabeled.tif'
files['MON2B'] = 'Hash_Segm_of_MON2B_20180723_ROI_Filtered_aniso_PVE_75_Tol_0.95_NumVS_1relabeled.tif'
files['MON3A'] = 'Hash_Segm_of_MON3A_20180723_ROI_Filtered_aniso_PVE_75_Tol_0.95_NumVS_1relabeled.tif'
files['MON3B'] = 'Hash_Segm_of_MON3B_20180723_ROI_Filtered_aniso_PVE_75_Tol_0.95_NumVS_1relabeled.tif'
path_save = '/home/german.martinez-carvajal/Bureau/These/Total_volumes_&_Global_fractions/manip_juillet'
chdir(path_save)
slices_start_deposit_layer = [157,693,180,660,292,106]
slices_end_deposit_layer = [2123,3013,3051,3228,3129,2889]
# Computing global fractions
for i in range(6):
sample = samples[i]
start = slices_start_deposit_layer[i]
end = slices_end_deposit_layer[i]
chdir(paths[sample])
print('Reading...')
image = tifffile.imread(files[sample])
compute_global_fractions(image, sample, start, end, len(image), path_save)
# # Computing Total Volumes
# for sample in samples:
# break
# chdir(paths[sample])
# image = tifffile.imread(files[sample])
# labels = [1,0,128,255] # 1 = outside the cylinder, 0 = voids, 128 = fouling material, 255 = gravel
# volumes = compute_volumes(image, labels, resolution) #
# volumes = np.append(volumes, np.sum(volumes[1:])) # do not consider outside voxels
# h,d,d = np.array(image.shape)*resolution/10 # in cm
# print(h,d,d)
# volume_cylinder = 0.25*np.pi*h*(d**2) #in cm3
# volumes = np.append(volumes, volume_cylinder)
# df[sample] = volumes
# chdir(path_save)
# resolution = 0.035 # mm/vox
# df = pd.DataFrame()
# df['Variable'] = ['Volume (cm3) of phase {}'.format(phase) for phase in ['outside', 'Voids', 'FM', 'Gravel', 'total', 'cylinder' ]]
# df.to_csv('Total_volumes.csv', sep = ';')
# print (df)
main()
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