# -*- coding: utf-8 -*- """ Created on Thu Mar 29 12:16:17 2018 @author: That guy that turns his pen all the fucking time. The purpose of this code is to make the spatial calibration for DaVis images. It consist in taking a doble-image, separate them in two and superimpose them. It also calculates the magnification in the picture [pixels]<->[mm] """ "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 datetime import cv2 import numpy as np import ReadIM import matplotlib.pyplot as plt import glob import os from skimage import exposure from scipy.optimize import curve_fit from functions_utile import draw_rectangle, image_correlation, test_homography, magnification """ Variables you can change in fonction of the experiment you are doing (f.e. the image to do the calibration with) """ working_directory = "C:\\Users\\Armando\\Desktop\\python_processing\\total_processing\\" nbr_of_experiments = 1 """ Variables initiation (don't touch unless you know what you are doing) Normally, this variables are the ones we will save to process the images later :) """ rectangle = {} rectangle['left'] = [] rectangle['right'] = [] """ This are the variables we need to process the spatial calibration and pH calibration, but once it is done we can delete them to save memory and hence go faster and farther !! """ images = {} file_name = [] pH_images ={} courbe = [] rectangle_final = [] pH_max_value_image = [] experiments = [] #Unlock this code if you want to process just one image, be wise though cause farther in the code #you have the sections to treat spatial and pH calibration #If you are working with DaViS here is where you write the name of the file #vbuff, vatts = ReadIM.extra.get_Buffer_andAttributeList(working_directory + 'calibration\\calibration_spatial\\' + 'mire.im7') #v_array, vbuff = ReadIM.extra.buffer_as_array(vbuff) #img_rescaled = exposure.rescale_intensity(v_array) #img_rescaled = img_rescaled[0] #images['image_mire'] = img_rescaled #If you are working with jpg,png or other more normal type of images #src_img = cv2.imread("E:\\Donnees\\ArmandoBulle\\python_treatement\\PLIF_test\\calibration\\calibration_spatial\\book1.jpg") #src_img_grey = cv2.cvtColor(src_img, cv2.COLOR_BGR2GRAY) #cv2.imshow('Grey Image',src_img_grey) #Test to see if the image you are using is the right one and is converted to grey #%% """This section of the code looks for all the files named ***pH***.im7 and ***mire***.im7 and stores them in image dictionary either with their pHname file eith with img_mire file. """ file_name = glob.glob(working_directory + "calibration\\**\\*.im7") img_rescaled = [] for e in range(len(file_name)): if file_name[e].find('pH') != -1: vbuff, vatts = ReadIM.extra.get_Buffer_andAttributeList(file_name[e]) v_array, vbuff = ReadIM.extra.buffer_as_array(vbuff) img_rescaled = exposure.rescale_intensity(v_array) img_rescaled = img_rescaled[0] name = file_name[e].rsplit('\\', 1)[1] images[name] = img_rescaled elif file_name[e].find('mire') != -1: vbuff, vatts = ReadIM.extra.get_Buffer_andAttributeList(file_name[e]) v_array, vbuff = ReadIM.extra.buffer_as_array(vbuff) img_rescaled = exposure.rescale_intensity(v_array) img_rescaled = img_rescaled[0] images['image_mire'] = img_rescaled #Cleaning the already used variables/ Memory clean up ReadIM.DestroyBuffer(vbuff) ReadIM.DestroyAttributeListSafe(vatts) del(img_rescaled, file_name ) #%% """ This section is the center of the spatial calibration, it takes a image, ask for the spatial points and creates the homography. It also crops the images and deduces the magnification (mm<->pixels). Finally it gives the superposed images. """ #Make two images out of the double-image #Obtain the image zone left images['image_left_mire'] = images['image_mire'].copy() images['image_left_mire'], rectangle['left'] = draw_rectangle(images['image_left_mire']) #Obtain the image zone right images['image_right_mire'] = images['image_mire'].copy() images['image_right_mire'], rectangle['right'] = draw_rectangle(images['image_right_mire']) #Get the magnification of the image magnification = magnification(images['image_left_mire']) #Obtain the homography between both images h , status = image_correlation(images['image_left_mire'],images['image_right_mire']) #Apply the homography transformation images['img_final'] = test_homography(h, status, images['image_left_mire'], images['image_right_mire']) cv2.namedWindow("Final Source Image",cv2.WINDOW_NORMAL) cv2.imshow("Final Source Image", images['img_final']) cv2.waitKey(0) cv2.destroyAllWindows() Super_position_images = images['img_final'].copy() #Cleaning the already used variables/ Memory clean up del(images['image_left_mire'],images['image_right_mire'],images['image_mire'], images['img_final']) #%% """pH image processing This is where the pH_images get normalized by the max value, which is get from the higher pH value The first boucle uses the spatial calibration to superpose the pH images. It also gets searches for the highest pH and hence the highest intensity. The second part divides the images by the highest pH image, to get values from 0 pH_max_value: pH_max_value = pH_value pH_max_value_image = pH_ratio del(pH_images_left, pH_images_right, pH, images, pH_ratio, titles) rectangle_final = draw_rectangle(pH_max_value_image)[1] cv2.waitKey(0) cv2.destroyAllWindows() pH_max_value_image = np.reciprocal(pH_max_value_image) for e in pH_images.keys(): pH_images[e] = np.divide(np.reciprocal(pH_images.get(e)),pH_max_value_image) pH_images[e] = pH_images[e][rectangle_final[0][1]:rectangle_final[1][1],rectangle_final[0][0]:rectangle_final[1][0]] if np.nanmean(pH_images[e], dtype = np.float64) <= 1: courbe.append([float(e),np.nanmean(pH_images[e], dtype = np.float64)]) elif np.nanmean(pH_images[e], dtype = np.float64) >1 and np.nanmean(pH_images[e], dtype = np.float64) <10 : courbe.append([float(e),1]) del(pH_images, pH_max_value, pH_max_value_image,) #%% courbe = np.array(courbe, dtype=np.float64) #courbe1 = np.array(courbe1, dtype=np.float64) courbe = np.sort(courbe,0) x1 = courbe[:,0] y1 = courbe[:,1] def tanh_fit(x, a, b, c, d): return a*(np.tanh(b*x+c)) + d def inverse_tanhfit(y ,a ,b ,c ,d): return (1/2*np.log((y+a-d)/(-y+a+d))-c)/(b) popt, pcov = curve_fit(tanh_fit, x1, y1,bounds=([0,0,-15,0], [1., 2.,5,1.])) fig = plt.figure() plt.title('Normalized intensity in function of pH') plt.xlabel('pH', color='plum') plt.ylabel('Normalized intensity (I/I\u2080)', color='0.5') # grayscale color plt.grid(b=True, which='major', color='k', linestyle='-') plt.grid(b=True, which='minor', color='k', linestyle='--') plt.plot(x1,y1, 'r.') x_full = np.arange(-1,15,0.01) y_full = np.arange(0,1,0.01) plt.plot(x_full, tanh_fit(x_full, *popt), '--',color = 'black',label='fit: A=%.3f, B=%.3f, C=%.3f D=%.3f \n A*tanh(B*pH+C) + D' % tuple(popt)) #plt.plot(y_full, inverse_tanhfit(y_full, *popt), '--',color = 'blue') plt.legend() plt.show() x_min = 0 x_max = courbe[len(courbe)-1][0] y_min = tanh_fit(x_min, *popt) y_max = tanh_fit(x_max, *popt) del(x1,y1,x_full) #%% #%% with open(working_directory + 'workfile.txt', mode='a+') as file: file.write('Printing the date recorded at %s.\n' % (datetime.datetime.now())) file.write('The homography matrix(3x3) is \n') [[file.write('%f \n' % x) for x in row]for row in h] file.write("This is the magnification %f \n" % magnification) file.write("This are the coordinates for image 1 crop upper corner : x = %d and y = %d \n" %(rectangle['left'][0][1],rectangle['left'][0][0])) file.write("lower corner : x = %d and y = %d \n" %(rectangle['left'][1][1],rectangle['left'][1][0])) file.write("This are the coordinates for image 2 crop upper corner : x = %d and y = %d \n" %(rectangle['right'][0][1],rectangle['right'][0][0])) file.write("lower corner : x = %d and y = %d \n" %(rectangle['right'][1][1],rectangle['right'][1][0])) file.write('The fitting curve with form : A*tanh(B*pH+C) + D \n') [[file.write('%f \n' % x) for x in popt]] file.write('------------------------------------*-------------------*------------------------------------\n') cv2.imwrite(working_directory + 'Super_position_images.tiff', Super_position_images) plt.savefig( working_directory +'Normalized intensity in function of pH.pdf') plt.savefig( working_directory + 'Normalized intensity in function of pH.png' ) #%% """ Traitement des données """ for root, dirs, files in os.walk(working_directory + "donnees"): for name in dirs: experiments.append(name) experiments_files = {} for e in range(len(experiments)): for root, dirs, files in os.walk(working_directory + "donnees\\"+ experiments[e]+ '\\'): experiments_files[experiments[e]] = files for i in range(len(experiments_files[experiments[e]])-1): vbuff, vatts = ReadIM.extra.get_Buffer_andAttributeList(working_directory + "donnees\\"+experiments[e]+'\\'+ experiments_files[experiments[e]][i]) v_array, vbuff = ReadIM.extra.buffer_as_array(vbuff) img_rescaled = exposure.rescale_intensity(v_array) image_experiment = img_rescaled[0] experiment_left = image_experiment.copy() experiment_left = experiment_left[rectangle['left'][0][1]:rectangle['left'][1][1],rectangle['left'][0][0]:rectangle['left'][1][0]] experiment_right = image_experiment.copy() experiment_right = experiment_right[rectangle['right'][0][1]:rectangle['right'][1][1],rectangle['right'][0][0]:rectangle['right'][1][0]] experiment_right = cv2.warpPerspective(experiment_right, h, (experiment_left.shape[1],experiment_left.shape[0])) experiment_ratio = np.divide(experiment_right,experiment_left) #experiment_ratio2 = np.divide(experiment_left,experiment_right) for y in range(len(experiment_ratio)-1): for x in range(len(experiment_ratio[0])-1): if np.isinf(experiment_ratio[y][x]): experiment_ratio[y][x] = 0 elif np.isnan(experiment_ratio[y][x]): experiment_ratio[y][x] = 0 elif experiment_ratio[y][x] < 1: experiment_ratio[y][x] = 1 #for y in range(len(experiment_ratio)-1): # for x in range(len(experiment_ratio[0])-1): # if np.isinf(experiment_ratio[y][x]): # experiment_ratio2[y][x] = 0 # elif np.isnan(experiment_ratio[y][x]): # experiment_ratio2[y][x] = 0 # elif experiment_ratio[y][x] < 1: # experiment_ratio2[y][x] = 1 path = working_directory + "donnees\\"+ 'results\\' + experiments[e] foto_name = path + '\\' + experiments_files[experiments[e]][i].rsplit('.')[0] #foto_name2 = path + '\\' + experiments_files[experiments[e]][i] + '2.tiff' if not os.path.exists(path): os.makedirs(path) #experiment_ratio = experiment_ratio/255. cv2.imshow('Grey Image1',experiment_ratio) experiment_ratio2 =np.reciprocal(experiment_ratio) cv2.imshow('Grey Image3',experiment_ratio) experiment_ratio=experiment_ratio.astype(np.float32) experiment_ratio=cv2.cvtColor(experiment_ratio , cv2.COLOR_GRAY2BGR) #cv2.imwrite(foto_name+ '.tiff', experiment_ratio) #cv2.imwrite(foto_name+ '.png', experiment_ratio) #cv2.imwrite(foto_name+ '.jpg', experiment_ratio) cv2.imshow('Grey Image2',experiment_ratio) #experiment_ratio2=experiment_ratio2.astype(np.float32) #experiment_ratio2=cv2.cvtColor(experiment_ratio2 , cv2.COLOR_GRAY2BGR) #cv2.imwrite(foto_name2, experiment_ratio)