# -*- 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