scientific_comp_projects/CODE/[python]thesis_old_scripts/image_processing_fit.py

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