324 lines
13 KiB
Python
324 lines
13 KiB
Python
# -*- 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<x<=1
|
|
"""
|
|
titles = list(images.items())
|
|
pH_max_value = 0
|
|
for e in range(len(titles)):
|
|
if titles[e][0].find('pH') != -1:
|
|
pH = titles[e][0].replace('.im7','').rsplit('pH', 1)[1]
|
|
pH_images[pH] = titles[e][1]
|
|
pH_images_left = pH_images[pH][rectangle['left'][0][1]:rectangle['left'][1][1],rectangle['left'][0][0]:rectangle['left'][1][0]]
|
|
pH_images_right = pH_images[pH][rectangle['right'][0][1]:rectangle['right'][1][1],rectangle['right'][0][0]:rectangle['right'][1][0]]
|
|
pH_images_right = cv2.warpPerspective(pH_images_right, h, (pH_images_left.shape[1],pH_images_left.shape[0]))
|
|
pH_ratio = np.divide(pH_images_right,pH_images_left)
|
|
for y in range(len(pH_ratio)-1):
|
|
for x in range(len(pH_ratio[0])-1):
|
|
if np.isinf(pH_ratio[y][x]):
|
|
pH_ratio[y][x] = 0
|
|
elif np.isnan(pH_ratio[y][x]):
|
|
pH_ratio[y][x] = 0
|
|
|
|
pH_images[pH] = pH_ratio
|
|
pH_value = float(pH)
|
|
if pH_value > 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)
|
|
|