# -*- coding: utf-8 -*-
"""
Created on Fri Jul 24 13:02:25 2020
@author: Sasha Marie Bakker
"""
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Choose the video you are tracking
video = "A00_01_G1"
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Import libraries
import numpy as np
import matplotlib.pyplot as plt
from skimage import io
import os
import statistics as stat
import pandas as pd
from numpy import diff, convolve
from scipy.stats import norm
def norm_dist_d(mu, sd):
"""Create normal distribition for convoluation."""
x = np.arange(-60,60,1)
npdf = norm.pdf(x,mu,sd)
dnpdf = diff(npdf)
return dnpdf
dnpdf = norm_dist_d(0,0.05)
pad = 0
"""
define the image dimensions and names
"""
if video == "A00_01_G1":
#a0, b0 = 0, 990
a0, b0 = 100, 855
y0, h0 = 35, 91
n_cols, n_rows = 402, 152
the_callibration_factor = 457.5750
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
elif video == "A01_01_G1":
#a0, b0 = 0, 805
a0, b0 = 190, 620
y0, h0 = 29, 86
n_cols, n_rows = 279, 142
the_callibration_factor = 462.1875
elif video == "A01_02_G1":
#a0, b0 = 0, 766
a0, b0 = 250, 685
y0, h0 = 24, 94
n_cols, n_rows = 257, 138
the_callibration_factor = 461.3625
elif video == "A01_03_G2":
#a0, b0 = 0, 792
a0, b0 = 260, 710
y0, h0 = 22, 100
n_cols, n_rows = 252, 146
the_callibration_factor = 462.1125
elif video == "A01_04_G2":
#a0, b0 = 0, 819
a0, b0 = 265, 655
y0, h0 = 26, 90
n_cols, n_rows = 256, 141
the_callibration_factor = 461.9125
elif video == "A01_05_G2":
#a0, b0 = 0, 774
a0, b0 = 155, 640
y0, h0 = 16, 105
n_cols, n_rows = 250, 134
the_callibration_factor = 462.3125
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
elif video == "A02_01_G3":
#a0, b0 = 0, 648
a0, b0 = 190, 525
y0, h0 = 12, 103
n_cols, n_rows = 222, 129
the_callibration_factor = 460.8125
elif video == "A02_02_G3":
#a0, b0 = 0, 560
a0, b0 = 130, 470
y0, h0 = 19, 99
n_cols, n_rows = 224, 128
the_callibration_factor = 461.2875
elif video == "A02_03_G3":
#a0, b0 = 0, 509
a0, b0 = 110, 445
y0, h0 = 17, 97
n_cols, n_rows = 221, 134
the_callibration_factor = 462.5000
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
elif video == "A03_01_G3":
#a0, b0 = 0, 513
a0, b0 = 85, 440
y0, h0 = 11, 106
n_cols, n_rows = 209, 126
the_callibration_factor = 459.6625
elif video == "A03_02_G3":
#a0, b0 = 0, 542
a0, b0 = 150, 445
y0, h0 = 16, 99
n_cols, n_rows = 203, 128
the_callibration_factor = 458.8520
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
elif video == "A04_01_G3":
#a0, b0 = 0, 338
a0, b0 = 95, 285
y0, h0 = 11, 105
n_cols, n_rows = 171, 126
the_callibration_factor = 462.0875
elif video == "A04_02_G3":
#a0, b0 = 0, 365
a0, b0 = 95, 300
y0, h0 = 15, 98
n_cols, n_rows = 173, 126
the_callibration_factor = 460.2750
elif video == "A04_03_G4":
#a0, b0 = 0, 341
a0, b0 = 95, 310
y0, h0 = 9, 103
n_cols, n_rows = 168, 126
the_callibration_factor = 459.0875
elif video == "A04_04_G4":
#a0, b0 = 0, 400
a0, b0 = 140, 350
y0, h0 = 17, 97
n_cols, n_rows = 178, 124
the_callibration_factor = 462.6500
elif video == "A04_05_G4":
#a0, b0 = 0, 329
a0, b0 = 75, 270
y0, h0 = 18, 96
n_cols, n_rows = 163, 133
the_callibration_factor = 462.2000
elif video == "A04_06_G4":
#a0, b0 = 0, 342
a0, b0 = 85, 310
y0, h0 = 14, 103
n_cols, n_rows = 179, 126
the_callibration_factor = 461.8000
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Nums = [str(item).zfill(4) for item in list(np.arange(a0,b0,1))]
n_imgs = len(Nums) # number of images
xs = np.arange(0,n_cols,1)
"""
define import and export paths
"""
intro_path = "C:/Users/sasha/Desktop/New_Test_Length_Series/Image_Series/"
load_image_path = intro_path + video
export_intensity_txt_path = intro_path + "Exports/" + video + "_Intensities/txt"
export_intensity_plots_path = intro_path + "Exports/" + video + "_Intensities/img"
#export_leading_minimum_plots_path = intro_path + "Exports/" + video + "_Intensities/cline"
#export_contact_minimum_plots_path = intro_path + "Exports/" + video + "_Intensities/ledge"
def export_intensity_plots():
j = 0
avgs_inten = np.zeros((n_imgs,n_cols))
sds_inten = np.zeros((n_imgs,n_cols))
avgs_convd = np.zeros((n_imgs,n_cols))
sds_convd = np.zeros((n_imgs,n_cols))
for N in Nums:
# condition puts this at 12 fps
if (int(N)%1)== 0:
os.chdir(load_image_path)
image_load_name = video + "_" + N + ".jpg"
image = io.imread(image_load_name)
image_red = image[:,:,0]
#print(image_red.shape)
avgs = list()
sds = list()
for col in range(len(image_red[0,:])):
y1, y2 = int(y0), int(y0 + h0 - 1)
avgs.append(stat.mean(image_red[y1:y2,col]))
sds.append(np.std(image_red[y1:y2,col]))
convd = convolve(avgs,dnpdf,'same')
convd_sds = convolve(sds, dnpdf, 'same')
# Uncomment next lines to create and save plots
"""
os.chdir(export_intensity_plots_path + "/inten")
plt.figure()
plt.title(image_load_name)
plt.xlabel("Columns")
plt.ylabel("Intensities")
plt.ylim(-200,300)
plt.grid()
plt.errorbar(xs,avgs,yerr=sds,ecolor='red')
plt.errorbar(xs,convd,yerr=convd_sds,ecolor='yellow')
plt.savefig("Intensities_" + video + "_" + N + ".jpg")
plt.close()
"""
#print(len(avgs_inten[j,:]),len(avgs))
avgs_inten[j,:] = np.array(avgs)
sds_inten[j,:] = np.array(sds)
avgs_convd[j,:] = np.array(convd)
sds_convd[j,:] = np.array(convd_sds)
j += 1
# Uncomment next lines to save data
#"""
os.chdir(export_intensity_txt_path)
np.savetxt("Average_Intensities_" + video + ".txt",avgs_inten)
np.savetxt("Stdev_Intensities_" + video + ".txt",sds_inten)
np.savetxt("Average_Derivatives_" + video + ".txt",avgs_convd)
np.savetxt("Stdev_Derivatives_" + video + ".txt",sds_convd)
#"""
print('done')
def find_leading_edge():
# Import text files
os.chdir(export_intensity_txt_path)
loaded = pd.read_csv("Average_Intensities_" + video + ".txt", sep=" ", header=None)
avgs_inten = loaded.to_numpy()
#loaded = pd.read_csv("Stdev_Intensities_" + video + ".txt", sep=" ", header=None)
#sds_inten = loaded.to_numpy()
loaded = pd.read_csv("Average_Derivatives_" + video + ".txt", sep=" ", header=None)
avgs_convd = loaded.to_numpy()
#loaded = pd.read_csv("Stdev_Derivatives_" + video + ".txt", sep=" ", header=None)
#sds_convd = loaded.to_numpy()
save_leading_edge = list()
for i in range(len(Nums)):
#for i in range(90,91):
N = int(Nums[i])
if (N%1) ==0:
derivative = avgs_convd[i,:]
intensity = avgs_inten[i,:]
magnify = derivative*xs
sort_magnify_indicies = np.argsort(magnify)
minimum_index = int(sort_magnify_indicies[0])
#reversed_derivative = derivative[::-1]
threshold = 10
for j in reversed(range(0,minimum_index)):
d_value = derivative[j]
if d_value > -threshold:
break
save_leading_edge.append(j)
# uncomment the next lines to create and save plots
"""
os.chdir(export_leading_minimum_plots_path)
plt.figure()
plt.grid()
plt.ylim(-300,300)
plt.title(f"Image {Nums[i]} from video {video}")
plt.plot(xs,intensity,color='magenta')
plt.plot(xs,derivative,color='orange')
#plt.plot(derivative*xs**2)
#plt.plot(minimum_index,derivative[minimum_index],'go')
plt.plot(j,d_value,'kx')
plt.plot(j,intensity[j],'kx')
plt.savefig("Minimums_" + video + "_" + str(N) + ".jpg")
plt.close()
#plt.show()
"""
os.chdir(export_intensity_txt_path)
np.savetxt("Leading_Edge_Column_Index_" + video + ".txt",save_leading_edge)
print('done')
def find_contact_line():
# Import text files
os.chdir(export_intensity_txt_path)
loaded = pd.read_csv("Leading_Edge_Column_Index_" + video + ".txt", sep=" ", header=None)
leading_edge = loaded.to_numpy()
loaded = pd.read_csv("Average_Intensities_" + video + ".txt", sep=" ", header=None)
avgs_inten = loaded.to_numpy()
loaded = pd.read_csv("Average_Derivatives_" + video + ".txt", sep=" ", header=None)
avgs_convd = loaded.to_numpy()
save_contact_line = list()
incorrect = list()
for i in range(len(Nums)):
#for i in range(len(Nums)-1,len(Nums)):
#for i in range(0,1):
N = int(Nums[i])
if (N%1) ==0:
derivative = avgs_convd[i,:]
intensity = avgs_inten[i,:]
leading_edge_index = int(leading_edge[i])
threshold = 25 # used 25
for j in range(len(derivative)):
value= derivative[j]
if value >= -threshold and value <= threshold:
derivative[j] = 0
if j >= leading_edge_index:
derivative[j] = 0
for j in reversed(range(0, leading_edge_index)):
value = derivative[j]
if value < 0:
break
k=4
hug = 0
mm=0
for q in reversed(range(0,j-k)):
mm+=1
hug = q
value = intensity[q]
if value > intensity[q+k]:
continue
else:
q=q+int(k/2)
#save_contact_line.append(q)
break
"""
# alternative
threshold = 5
for q in reversed(range(0,j-k)):
mm+=1
hug = q
value = derivative[q]
if value > -threshold:
break
"""
#max_col = np.argsort(intensity[q:q+k])[::-1][0]+q
save_contact_line.append(hug) # q otherwise
if mm == 0:
incorrect.append(i)
os.chdir(export_intensity_txt_path)
np.savetxt("Contact_Line_Column_Index_" + video + ".txt",save_contact_line)
#os.chdir(intro_path + "New_Intensity_Data/" + "003" + "_Intensities")
#loaded = pd.read_csv("Leading_Edge_Column_Index_" + "003" + ".txt", sep=" ", header=None)
leading_edge_003 = leading_edge
#loaded = pd.read_csv("Contact_Line_Column_Index_" + "003" + ".txt", sep=" ", header=None)
loaded = pd.read_csv("Contact_Line_Column_Index_" + video + ".txt", sep=" ", header=None)
contact_line_003 = loaded.to_numpy()
incorrect = np.array(incorrect)
if len(incorrect) > 0:
contact_line_003 = np.delete(contact_line_003,incorrect)
leading_edge_003 = np.delete(leading_edge_003,incorrect)
callibrate = the_callibration_factor
diff_003 = leading_edge_003 - contact_line_003
diff_003 = diff_003/callibrate
plt.figure()
plt.ylim(0,2.5)
x = np.arange(0,len(leading_edge_003),1,)/60
plt.plot(x,leading_edge_003/callibrate,color='green',label='leading edge')
plt.plot(x,contact_line_003/callibrate,color='magenta',label='contact line')
plt.title(f"video {video}")
plt.xlabel("Time (s)")
plt.ylabel("Distance (cm)")
plt.grid()
plt.legend()
plt.show()
x = np.arange(0,len(diff_003),1,)/60
"""
wlist = list()
for w in range(0,len(diff_003)):
val = diff_003[w]
if val <= 0.06:
wlist.append(w)
warr = np.array(wlist)
diff_003 = np.delete(diff_003,warr)
x = np.delete(x,warr)
print(warr)
"""
plt.figure()
plt.plot(x, diff_003,'kx')
m,b = np.polyfit(x[:],diff_003[:],1)
plt.plot(x[:],m*x[:]+b,color='red')
plt.title(f"video {video}, m = {m}, b = {b}")
plt.xlabel("Time (s)")
plt.ylabel("Distance (cm)")
plt.grid()
plt.show()