517 lines
22 KiB
Python
517 lines
22 KiB
Python
# -*- coding: utf-8 -*-
|
|
"""
|
|
This script is made to process images for Quenched-PLIF method in fluid dynamics.
|
|
Before using it, you must have already obtained the fit values (a,b,c and d) for the inverse hyperbolic tangent.
|
|
Also, the buffers (images from DaVis* Lavision) must already be spatially calibrated and in ratio (divided)
|
|
|
|
Finally, the script is divided into three major sections :
|
|
-the first one is the initialization, this section contains the libraries to be imported for the script,
|
|
the defined functions (such as the fitting courbe or the animation and saving into images) and the working directory definition
|
|
|
|
-the second one is the treatment in itself, this part of the code is made to treat all the buffers that appear in the
|
|
experiments list.
|
|
|
|
-the third one is the test part of the script, this is somewhat like the second part but rather than processing all the
|
|
buffers, you have to define just one image to process. This helps to see the outcome and modify the second script before running it
|
|
|
|
|
|
Created on Mon May 14 18:14:28 2018
|
|
|
|
@author: Armando FEMAT ORTIZ, PhD Student at Ecole Centrale Lyon, Fluid Mechanics and Acoustics Laboratory (LMFA)
|
|
"""
|
|
|
|
"We tried to do this as cleanly as possible so that some PhD student in the future"
|
|
"won't get on a plane at 2:00am and come hunt us down with a crowbar."
|
|
|
|
|
|
import os
|
|
import ReadIM
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
import matplotlib.animation as animation
|
|
import glob
|
|
from matplotlib.colors import LogNorm
|
|
from matplotlib.animation import FFMpegWriter
|
|
#from skimage.filters.rank import mean_bilateral
|
|
#from skimage.morphology import disk
|
|
#import scipy.misc
|
|
|
|
#Courbe fit values - as well as the maximum and minimum values for the inverse fit function
|
|
global a, b, c, d
|
|
a = 0.316617
|
|
b = 1.26288
|
|
c = -7.70744
|
|
d = 0.6722
|
|
|
|
|
|
y_min = 0.3559
|
|
y_max = 0.98878
|
|
|
|
|
|
|
|
|
|
|
|
#Numerical to physical equivalences
|
|
nappe_laser = 0.250 #Thickness of laser en mm
|
|
dimension_pixel = 0.026244 #Equivalence mm/px
|
|
|
|
#Conversion to SI units
|
|
dimension_pixel_m = dimension_pixel * 10**(-3) #Equivalence mm/px
|
|
nappe_laser_m = nappe_laser * 10**(-3) #Thickness of laser en m
|
|
surface_pixel = dimension_pixel_m**2 #Surface equivalence of a pixel (m^2)
|
|
|
|
#The working directory, where all the images are stored
|
|
working_directory = 'D:\\python_processing\\imagesV4\\'
|
|
plt.ioff()
|
|
#help to accelerate script by not showing images
|
|
#Plotting options
|
|
global interpolation, vmax_value, vmin_value, colormap
|
|
colormap = plt.get_cmap('viridis')
|
|
vmin_value = 0.00001 #Minimum value for the pixel value of a processed image
|
|
vmax_value = 40 #Maximum value for the pixel value of a processed image
|
|
interpolation = 'bicubic' #Output image interpolation (see matplotlib interpolations for more information)
|
|
|
|
#Fit-inverse function
|
|
def inverse_tanhfit(y):
|
|
return (1/2*np.log((a+y-d)/(a-y+d))-c)/(b)
|
|
|
|
#Animation function
|
|
def animate(i):
|
|
arr = frames[i]
|
|
im.set_interpolation(interpolation)
|
|
im.set_data(arr)
|
|
vmin_value
|
|
vmax_value
|
|
im.set_clim(vmin_value, vmax_value)
|
|
im.set_cmap(colormap)
|
|
tx.set_text('Frame {0}'.format(i))
|
|
plt.savefig(newpath + 'Frame {0}'.format(i) + ".png",bbox_inches='tight')
|
|
# In this version you don't have to do anything to the colorbar,
|
|
# it updates itself when the mappable it watches (im) changes
|
|
|
|
experiments = []
|
|
for root, dirs, files in os.walk(working_directory):
|
|
for name in dirs:
|
|
experiments.append(name)
|
|
break #exit after first subdirectories
|
|
|
|
del(dirs,files,root,name)
|
|
|
|
#%%In this section only the files are treated, no images are created hence the treatment time is much lower
|
|
|
|
#Data to be modified
|
|
file_type = '.im7' #The file type used (if it is not a DaVis file-type the code must be modified)
|
|
cut_image_y0 = 0 #The borders that will be cropped from the image [y0,yf,x0,xf]
|
|
cut_image_yf = 0 #xf and yf, are the values to be rested from the total images(hence if 0 the max values are used)
|
|
cut_image_x0 = 0
|
|
cut_image_xf = 0
|
|
radius = 0
|
|
r = []
|
|
radii = False
|
|
boundary_two = 0
|
|
boundary_one = 0
|
|
|
|
|
|
#max_value = []
|
|
elements_number_max1 = []
|
|
elements_number_max2 = []
|
|
#min_value = []
|
|
#elements_number_min = []
|
|
massCO2_numer_max = {}
|
|
|
|
|
|
for e in range(len(experiments)):
|
|
|
|
#Search for all the images (they should finish in *.im7, if not change the) in the 'e' experiment
|
|
file_name = glob.glob(working_directory +experiments[e]+ "\\*"+file_type)
|
|
mass_co2 = []
|
|
newpath = working_directory +experiments[e]+ "\\resultstest\\"
|
|
if not os.path.exists(newpath):
|
|
os.makedirs(newpath)
|
|
#Looping for all images in a given experiment
|
|
|
|
for i in range(len(file_name)):
|
|
|
|
#Reading and tranforming the buffer image of DaVis into a numpy array
|
|
vbuff, vatts = ReadIM.extra.get_Buffer_andAttributeList(file_name[i])
|
|
v_array, vbuff = ReadIM.extra.buffer_as_array(vbuff)
|
|
v_array = v_array[0][cut_image_y0:len(v_array[0])-cut_image_yf,cut_image_x0:len(v_array[0][0])-cut_image_xf]
|
|
modif_array = np.array(v_array,dtype = float)
|
|
original_array = np.array(v_array,dtype = float)
|
|
del(vbuff, v_array, vatts)
|
|
|
|
#Looping to apply the inverse-hyperbolic tangent function
|
|
#We get intensity as the input value, and pH as output value.
|
|
for y in range(len(original_array)):
|
|
|
|
begin = False
|
|
inside = False
|
|
end = False
|
|
|
|
for x in range(len(original_array[0])):
|
|
#If the value of the intensity is between ymax and ymin we obtain the pH through the 'inverse_tanhfit' function defined above
|
|
|
|
if (modif_array[y][x] == 0) and (begin) and (not end):
|
|
boundary_two = ( x - 1 )
|
|
inside = False
|
|
end = True
|
|
|
|
elif (modif_array[y][x] != 0) and (not inside):
|
|
boundary_one = x
|
|
begin = True
|
|
inside = True
|
|
|
|
if (original_array[y][x] <= y_max) and (original_array[y][x] >= y_min):
|
|
#The CO2 concentration is then obtained from the pH value
|
|
#modif_array[y][x] = 207608319.9386*np.power(inverse_tanhfit(original_array[y][x]),-9.7399) #Concentration CO2 en g/m^3 (equal to mg/l)(from measures made by Tom LACASSAGNE)
|
|
modif_array[y][x] = 115467070469*np.exp(-4.63*inverse_tanhfit(original_array[y][x]))*1.1 #Concentration CO2 en g/m^3 (equal to mg/l)(numerical values)
|
|
|
|
|
|
m = np.split(modif_array[y],[boundary_one,boundary_two])
|
|
zeros = len(modif_array[0])-len(m[1])
|
|
modif_array[y] = np.concatenate((np.zeros(zeros//2),m[1],np.zeros(int(np.ceil(zeros/2)))))
|
|
|
|
modif_array[modif_array>40] = 0
|
|
|
|
#The values under 9.45e-06 and over 100 (g/m^3) are considered noise and hence turned to 0
|
|
#modif_array[modif_array<9.45e-06] = 0
|
|
|
|
'''
|
|
for y in range(len(modif_array)):
|
|
|
|
begin = False
|
|
inside = False
|
|
end = False
|
|
|
|
for x in range(len(modif_array[0])):
|
|
if (modif_array[y][x] == 0) & (begin == True) & (end == False):
|
|
boundary_two = ( x - 1 )
|
|
inside = False
|
|
end = True
|
|
|
|
elif (modif_array[y][x] != 0) & (inside == False):
|
|
boundary_one = x
|
|
begin = True
|
|
inside = True
|
|
|
|
m = np.split(modif_array[y],[boundary_one,boundary_two])
|
|
zeros = len(modif_array[0])-len(m[1])
|
|
modif_array[y] = np.concatenate((np.zeros(zeros//2),m[1],np.zeros(int(np.ceil(zeros/2)))))
|
|
|
|
#Count the number of pixles in a range of values(dont forget to erase the # before the for loop)
|
|
max_value_actuel = np.amax(modif_array)
|
|
max_value.append(max_value_actuel)
|
|
elements_number_max.append(
|
|
len(np.where((modif_array >= (max_value_actuel - 2.0) ) & (modif_array <=max_value_actuel))[0]))
|
|
|
|
min_value_actuel = np.amin(modif_array)
|
|
min_value.append(min_value_actuel)
|
|
elements_number_min.append(
|
|
len(np.where((modif_array >= min_value_actuel) & (modif_array <=(min_value_actuel)))[0]))
|
|
'''
|
|
|
|
#Here we cut the images in two and inverse the second one.
|
|
cut_image1 = modif_array[0:len(modif_array),0:len(modif_array[0])//2]
|
|
cut_image2 = modif_array[0:len(modif_array),len(modif_array[0])//2:len(modif_array[0])]
|
|
cut_image2 = np.fliplr(cut_image2[0:len(cut_image1),0:len(cut_image1[0])])
|
|
|
|
#Then we take the mean value for both images.
|
|
modif_array = np.mean(np.array([cut_image1,cut_image2]), axis=0,dtype=np.float64)
|
|
#Error checking
|
|
modif_array = np.concatenate((modif_array,modif_array[0:len(modif_array),len(modif_array[0])-2:len(modif_array[0])]), axis = 1)
|
|
#modif_array = modif_array[:len(modif_array),:len(modif_array[0])-2]
|
|
|
|
|
|
#This section, creates the r-vector (which have the same size as cutted images)
|
|
if radii == False:
|
|
r = []
|
|
radius = dimension_pixel_m
|
|
radii = True
|
|
|
|
for k in range(len(modif_array[0])):
|
|
r.append(radius)
|
|
radius += dimension_pixel_m
|
|
|
|
r = np.flipud(r)
|
|
|
|
#Multiplication of the image [C] by the {r} radius vector
|
|
modif_array = (modif_array) @ r
|
|
|
|
#Final calculation to obtain the mass of CO2 in the image
|
|
mass_co2.append((surface_pixel * np.pi * 2 * modif_array.sum())/nappe_laser_m) #The mass is obtained in (g)
|
|
|
|
massCO2_numer_max[experiments[e]] = mass_co2
|
|
radii = False
|
|
|
|
#%% In this section all the folders with images which are referenced in the list 'experiments' will be treated.
|
|
# The variables before the start of the 'for' loop are to be adapted for each set of images, hence testes should be made first
|
|
|
|
#Data to be modified
|
|
metadata = dict(title='', artist='Armando FEMAT',
|
|
comment='Movie support!')
|
|
file_type = '.im7' #The file type used (if it is not a DaVis file-type the code must be modified)
|
|
frequency = 2 #This is the frequency of data adquisition (for the video to be in real-time)
|
|
cut_image_y0 = 0 #The borders that will be cropped from the image [y0,yf,x0,xf]
|
|
cut_image_yf = 0 #xf and yf, are the values to be rested from the total images(hence if 0 the max values are used)
|
|
cut_image_x0 = 0
|
|
cut_image_xf = 0
|
|
radius = 0
|
|
r=[]
|
|
mass_co2 = []
|
|
|
|
|
|
|
|
#Initialization variables
|
|
fig = plt.figure()
|
|
|
|
for e in range(len(experiments)):
|
|
|
|
#Search for all the images (they should finish in *.im7, if not change the) in the 'e' experiment
|
|
file_name = glob.glob(working_directory +experiments[e]+ "\\*"+file_type)
|
|
metadata['title'] = experiments[e]
|
|
writer = FFMpegWriter(fps=frequency, metadata=metadata)
|
|
|
|
#Creates the 'results' directory where all the processed images will be stored
|
|
newpath = working_directory +experiments[e]+ "\\resultstest\\"
|
|
if not os.path.exists(newpath):
|
|
os.makedirs(newpath)
|
|
|
|
#Restart and creation of the plotting for matplotlib
|
|
frames = []
|
|
ax = fig.add_subplot(111)
|
|
|
|
#Looping for all images in a given experiment
|
|
for i in range(len(file_name)):
|
|
#Reading and tranforming the buffer image of DaVis into a numpy array
|
|
vbuff, vatts = ReadIM.extra.get_Buffer_andAttributeList(file_name[i])
|
|
v_array, vbuff = ReadIM.extra.buffer_as_array(vbuff)
|
|
v_array = v_array[0][cut_image_y0:len(v_array[0])-cut_image_yf,cut_image_x0:len(v_array[0][0])-cut_image_xf]
|
|
modif_array = np.array(v_array,dtype = float)
|
|
original_array = np.array(v_array,dtype = float)
|
|
|
|
#Looping to apply the inverse-hyperbolic tangent function
|
|
#We get intensity as the input value, and pH as output value.
|
|
for y in range(len(original_array)-1):
|
|
for x in range(len(original_array[0])-1):
|
|
#If the value of the intensity is really big (f.e. due to reflexion) we impose a pH of 12
|
|
if original_array[y][x] >= y_max:
|
|
modif_array[y][x] = 9
|
|
#Similar, if the value is small, we impose a pH of 4
|
|
elif original_array[y][x] <= y_min:
|
|
modif_array[y][x] = 3.9
|
|
#Otherwise, we obtain the pH through the 'inverse_tanhfit' function defined above
|
|
else:
|
|
modif_array[y][x] = inverse_tanhfit(original_array[y][x])*0.9
|
|
|
|
#The CO2 concentration is then obtained from the pH value
|
|
#modif_array = 207608319.9386*np.power(modif_array,-9.7399) #Concentration CO2 en g/m^3 (equal to mg/l)(from measures made by Tom LACASSAGNE)
|
|
modif_array = 115467070469*np.exp(-4.63*modif_array) #Concentration CO2 en g/m^3 (equal to mg/l)(theoretical values)
|
|
|
|
#The values under 0.1 and over 40 (g/m^3) are considered noise and hence turned to 0
|
|
modif_array[modif_array<9.45e-06] = 0
|
|
modif_array[modif_array>10] = 0
|
|
|
|
#This line 'appends' the images in the list 'frames' to create the video afterwards
|
|
frames.append(modif_array)
|
|
|
|
#Here we cut the images in two and inverse the second one.
|
|
cut_image1 = modif_array[0:len(modif_array),0:len(modif_array[0])//2]
|
|
cut_image2 = np.fliplr(modif_array[0:len(modif_array),len(modif_array[0])//2:len(modif_array[0])-1])
|
|
|
|
#Then we take the mean value for both images.
|
|
modif_array = np.mean(np.array([cut_image1,cut_image2]), axis=0,dtype=np.float64)
|
|
|
|
#This section, creates the r-vector (which have the same size as cutted images)
|
|
if r == []:
|
|
for k in range(len(modif_array[0])):
|
|
radius += dimension_pixel_m
|
|
r.append(radius)
|
|
#Multiplication of the image [C] by the {r} radius vector
|
|
modif_array = (modif_array) @ r
|
|
|
|
#Final calculation to obtain the mass of CO2 in the image
|
|
mass_co2.append((surface_pixel * np.pi * 2 * modif_array.sum())/nappe_laser_m) #The mass is obtained in (g)
|
|
|
|
#Plot for the first image (helps to intialize the plotting function)
|
|
cv0 = frames[0]
|
|
im = ax.imshow(cv0,interpolation =interpolation , cmap=colormap , origin='upper', animated = True)
|
|
|
|
tx = ax.set_title('Frame 0')
|
|
plt.savefig(newpath + 'Frame {0}'.format(i) + ".png",bbox_inches='tight')
|
|
|
|
#Saving the images into *.png format and the video in a *.mp4 format
|
|
ani = animation.FuncAnimation(fig, animate, frames=len(frames),interval=1/frequency*1000, repeat = False)
|
|
ani.save(newpath + experiments[e]+".mp4", writer=writer)
|
|
|
|
fig.clf()
|
|
|
|
plt.close('all')
|
|
|
|
#%% This is the last section for image processing. Here is the laboratory to tune and upgrade the processing parameters.
|
|
#It is divided into two sections, the first one is to choose the image as well as the borders for the image (start with 0 here, twitch afterwards)
|
|
|
|
cut_image_y0 = 150 #The borders that will be cropped from the image [y0,yf,x0,xf]
|
|
cut_image_yf = 300 #xf and yf, are the values to be rested from the total images(hence if 0 the max values are used)
|
|
cut_image_x0 = 0
|
|
cut_image_xf = 0
|
|
radius = 0
|
|
r=[]
|
|
|
|
|
|
#The experiment and image to process (only one of each will be used since it is a test)
|
|
working_dir = 'D:\\python_processing\\imagesV2\\eau\\traite\\'
|
|
image_to_process = 'B00021.im7'
|
|
|
|
|
|
|
|
vbuff, vatts = ReadIM.extra.get_Buffer_andAttributeList(working_dir+ image_to_process)
|
|
v_array, vbuff = ReadIM.extra.buffer_as_array(vbuff)
|
|
original_array = v_array[0][cut_image_y0:len(v_array[0])-cut_image_yf,cut_image_x0:len(v_array[0][0])-cut_image_xf]
|
|
original_array = np.array(original_array,dtype = float)
|
|
modif_array = np.array(original_array,dtype = float)
|
|
for y in range(len(original_array)-1):
|
|
for x in range(len(original_array[0])-1):
|
|
if original_array[y][x] >= y_max:
|
|
modif_array[y][x] = 8
|
|
elif original_array[y][x] <= y_min:
|
|
modif_array[y][x] = 4
|
|
else:
|
|
modif_array[y][x] = inverse_tanhfit(original_array[y][x])
|
|
|
|
#The CO2 concentration is then obtained from the pH value
|
|
#modif_array = 207608319.9386*np.power(modif_array,-9.7399) #Concentration CO2 en g/m^3 (equal to mg/l)(from measures made by Tom LACASSAGNE)
|
|
modif_array = 115467070469*np.exp(-4.63*modif_array) #Concentration CO2 en g/m^3 (equal to mg/l)(theoretical values)
|
|
modif_array[modif_array<8e-13] = 0
|
|
modif_array[modif_array>1] = 0
|
|
|
|
#The values under 0.1 and over 40 (g/m^3) are considered noise and hence turned to 0
|
|
cut_image1 = modif_array[0:len(modif_array),0:len(modif_array[0])//2]
|
|
cut_image2 = np.fliplr(modif_array[0:len(modif_array),len(modif_array[0])//2:len(modif_array[0])-1])
|
|
|
|
#Then we take the mean value for both images.
|
|
modif_array2 = np.mean(np.array([cut_image1,cut_image2]), axis=0,dtype=np.float64)
|
|
|
|
#This section, creates the r-vector (which have the same size as cutted images)
|
|
if r == []:
|
|
for k in range(len(modif_array2[0])):
|
|
radius += dimension_pixel_m
|
|
r.append(radius)
|
|
|
|
#Multiplication of the image [C] by the {r} radius vector
|
|
modif_array2 = (modif_array2) @ r
|
|
|
|
#Final calculation to obtain the mass of CO2 in the image
|
|
mass_co2 = (surface_pixel * np.pi * 2 * modif_array2.sum())/nappe_laser_m #The mass is obtained in (g)
|
|
|
|
|
|
#%%This is the second section that will be used to plot the data obtained
|
|
|
|
#Use this configuration if you don't know the image exposure scale (the CO2 concentration in our case)
|
|
'''
|
|
lev_exp = np.arange(np.floor(np.log10(modif_array.min())-1), \
|
|
np.ceil(np.log10(modif_array.max())+1))
|
|
levs = np.power(10, lev_exp)
|
|
minima = levs[0]
|
|
maxima = levs[-1]
|
|
'''
|
|
#If you have an idea of the minima and maximum, please use this configuration and adjust the parameters
|
|
minima = 0.000007
|
|
maxima = 0.5
|
|
|
|
|
|
#Here is where the magic happens, thank you matplotlib!
|
|
fig = plt.figure()
|
|
a1 = fig.add_subplot(1, 2, 1)
|
|
|
|
#Plotting of the image
|
|
im = plt.imshow(modif_array, interpolation=interpolation, cmap=colormap,origin='upper',norm = LogNorm(minima , maxima))
|
|
plt.ylabel('y (mm)')
|
|
plt.xlabel('x (mm)')
|
|
plt.colorbar(im, orientation='horizontal')
|
|
#norm = colors.LogNorm(vmin=minima, vmax=maxima)
|
|
|
|
|
|
a2 = fig.add_subplot(1, 2, 2)
|
|
|
|
#Plotting of the image
|
|
#im = plt.imshow(original_array, interpolation=interpolation, cmap='gray',origin='upper',vmin=minima, vmax=maxima)
|
|
|
|
histo = plt.hist(modif_array, bins=10**np.linspace(np.log10(minima), np.log10(maxima), 10), fc='k', ec='k' )
|
|
|
|
#
|
|
#**np.linspace(np.log10(minima), np.log10(maxima), 10)
|
|
plt.gca().set_xscale('log')
|
|
a2.set_title('Compare and contrast')
|
|
|
|
|
|
plt.tight_layout()
|
|
plt.show()
|
|
#plt.savefig(working_dir + name + ".png",bbox_inches='tight')
|
|
|
|
|
|
#%%
|
|
|
|
|
|
#Data to be modified
|
|
file_type = '.im7' #The file type used (if it is not a DaVis file-type the code must be modified)
|
|
|
|
|
|
position_final = {}
|
|
|
|
|
|
for e in range(len(experiments)):
|
|
|
|
#Search for all the images (they should finish in *.im7, if not change the) in the 'e' experiment
|
|
file_name = glob.glob(working_directory +experiments[e]+ "\\*"+file_type)
|
|
position = []
|
|
|
|
for i in range(len(file_name)):
|
|
|
|
#Reading and tranforming the buffer image of DaVis into a numpy array
|
|
vbuff, vatts = ReadIM.extra.get_Buffer_andAttributeList(file_name[i])
|
|
v_array, vbuff = ReadIM.extra.buffer_as_array(vbuff)
|
|
modif_array = np.array(v_array[0],dtype = float)
|
|
original_array = np.array(v_array[0],dtype = float)
|
|
del(vbuff, v_array, vatts)
|
|
|
|
#Looping to apply the inverse-hyperbolic tangent function
|
|
#We get intensity as the input value, and pH as output value.
|
|
y_position = 0
|
|
done1 = False
|
|
count = 0
|
|
for y in range(len(original_array)):
|
|
|
|
begin = False
|
|
inside = False
|
|
end = False
|
|
|
|
for x in range(len(original_array[0])):
|
|
#If the value of the intensity is between ymax and ymin we obtain the pH through the 'inverse_tanhfit' function defined above
|
|
|
|
if (modif_array[y][x] == 0) & (begin == True) & (end == False):
|
|
boundary_two = ( x - 1 )
|
|
inside = False
|
|
end = True
|
|
|
|
elif (modif_array[y][x] != 0) & (inside == False):
|
|
boundary_one = x
|
|
begin = True
|
|
inside = True
|
|
done1 = True
|
|
|
|
if (done1 == True) & (count != 0) & (begin == True):
|
|
count = 0
|
|
|
|
if (done1 == True) & (begin == False) & (count == 0):
|
|
y_position = y
|
|
count += 1
|
|
|
|
if (done1 == True) & (begin == False) & (count != 0):
|
|
count += 1
|
|
|
|
if (done1 == True) & (begin == False) & count==10:
|
|
break
|
|
position.append(y_position) #The mass is obtained in (g)
|
|
|
|
position_final[experiments[e]] = position
|