168 lines
6.9 KiB
Python
Executable File
168 lines
6.9 KiB
Python
Executable File
##!/usr/bin/env python3
|
|
import numpy as np
|
|
import cv2
|
|
from pathlib import Path
|
|
import math
|
|
import h5py
|
|
import time
|
|
import scipy.stats
|
|
|
|
|
|
def create_slope_image(image_width: int, image_height: int, p1_angle: float, p1_radius: float, bg_color: float = 0.0, fg_color: float = 1.0):
|
|
"""
|
|
creates an image containing 2 regions separated by a line defined as the tangent of a point p1 on a centered circle
|
|
|
|
the semi-plane the contains the origin of the circle is colored with fg_color, while the other semi plane is assigned pixel values bg_color
|
|
"""
|
|
slope_image = np.ndarray(shape=(image_height, image_width), dtype=float)
|
|
# work out the equation of the tangent in the form a.x + b.y + c = 0
|
|
c = np.array([image_height * 0.5, image_width * 0.5])
|
|
p1 = c + np.array([math.cos(p1_angle), math.sin(p1_angle)]) * p1_radius
|
|
cp1 = p1 - c
|
|
a = cp1[0]
|
|
b = cp1[1]
|
|
c = - p1[0] * cp1[0] - p1[1] * cp1[1]
|
|
x = np.tile(np.arange(image_width).reshape((1, image_width)), (image_height, 1))
|
|
y = np.tile(np.arange(image_height).reshape((image_height, 1)), (1, image_width))
|
|
# print('x=', x)
|
|
# print('y=', y)
|
|
retval, slope_image = cv2.threshold(x * a + y * b + c, 0.0, fg_color, type=cv2.THRESH_BINARY_INV)
|
|
# print('slope_image=', slope_image)
|
|
cv2.imwrite('slope.tif', slope_image)
|
|
return slope_image
|
|
|
|
|
|
class TracData():
|
|
|
|
def __init__(self):
|
|
self.pts = None
|
|
self.th = None
|
|
self.num_frames = 0
|
|
self.image_size = None
|
|
|
|
|
|
def trac_data_to_hdf5(trac_data: TracData, hdf5_file_path: Path):
|
|
version = '3.0 (22/05/2021) | J. Heyman'
|
|
f = h5py.File(hdf5_file_path, 'w')
|
|
f.attrs['version'] = 'HDF5 file made with TracTrac Python v' + version
|
|
f.attrs['date'] = time.strftime("%d/%m/%Y")
|
|
f.attrs['nFrames'] = trac_data.num_frames
|
|
f.attrs['size'] = trac_data.image_size
|
|
for items in trac_data.th.keys():
|
|
f.attrs['th:' + items] = trac_data.th[items]
|
|
f.create_dataset("Frame", data=np.uint16(trac_data.pts[:, 0]))
|
|
f.create_dataset("Id", data=np.uint16(trac_data.pts[:, 1]))
|
|
f.create_dataset("x", data=np.float32(trac_data.pts[:, 2:4]))
|
|
f.create_dataset("u", data=np.float32(trac_data.pts[:, 4:6]))
|
|
f.create_dataset("a", data=np.float32(trac_data.pts[:, 6:8]))
|
|
f.create_dataset("u_motion", data=np.float32(trac_data.pts[:, 8:10]))
|
|
f.create_dataset("err_motion", data=np.uint32(trac_data.pts[:, 10]))
|
|
f.close()
|
|
|
|
|
|
def hdf5_to_trac_data(hdf5_file_path: Path):
|
|
|
|
trac_data = TracData()
|
|
f = h5py.File(hdf5_file_path, 'r')
|
|
|
|
version = f.attrs['version'] # noqa
|
|
date = f.attrs['date'] # noqa
|
|
trac_data.num_frames = f.attrs['nFrames']
|
|
trac_data.image_size = f.attrs['size']
|
|
print(type(trac_data.image_size))
|
|
# f.attrs['size'] = I0f.shape
|
|
# for items in th[0].keys():
|
|
# f.attrs['th:' + items] = th[0][items]
|
|
num_pts = f['Frame'].shape[0]
|
|
trac_data.pts = np.ndarray(shape=(num_pts, 11), dtype=float)
|
|
trac_data.pts[:, 0] = f['Frame']
|
|
trac_data.pts[:, 1] = f['Id']
|
|
trac_data.pts[:, 2:4] = f['x']
|
|
trac_data.pts[:, 4:6] = f['u']
|
|
trac_data.pts[:, 6:8] = f['a']
|
|
trac_data.pts[:, 8:10] = f['u_motion']
|
|
trac_data.pts[:, 10] = f['err_motion']
|
|
f.close()
|
|
return trac_data
|
|
|
|
|
|
class SlopeFinder():
|
|
|
|
def __init__(self, beads_radius: float):
|
|
self.beads_radius = beads_radius
|
|
|
|
def find_slope(self, trac_data: TracData, frame_index: int):
|
|
isbead_image = SlopeFinder.create_isbead_image(trac_data, frame_index, self.beads_radius)
|
|
cv2.imwrite('isbead_%04d.tif' % frame_index, isbead_image)
|
|
if False:
|
|
kernel = np.ones((15, 15), np.uint8)
|
|
closing = cv2.morphologyEx(isbead_image, cv2.MORPH_CLOSE, kernel)
|
|
cv2.imwrite('closing_%04d.tif' % frame_index, closing)
|
|
# remove the jumping beads
|
|
# apply connected component analysis to the thresholded image
|
|
# thresh = cv2.threshold(isbead_image, 0, 255, cv2.THRESH_BINARY_INV | cv2.THRESH_OTSU)[1]
|
|
thresh = cv2.compare(isbead_image, 0.5, cv2.CMP_GT)
|
|
connectivity = 4
|
|
output = cv2.connectedComponentsWithStats(thresh, connectivity, cv2.CV_32S)
|
|
(num_labels, labels_image, stats, centroids) = output
|
|
print('num_labels: ', num_labels)
|
|
print('labels_image: ', labels_image)
|
|
print('stats: ', stats)
|
|
print('centroids: ', centroids)
|
|
cv2.imwrite('labels_%04d.tif' % frame_index, labels_image)
|
|
is_non_jumping_bead_image = np.zeros(shape=labels_image.shape, dtype=np.uint8)
|
|
bead_area = math.pi * self.beads_radius * self.beads_radius
|
|
for label in range(1, num_labels): # ignore the first label, as it's the background
|
|
area = stats[label][4]
|
|
if area > bead_area * 1.1:
|
|
is_non_jumping_bead_image[labels_image == label] = 255
|
|
cv2.imwrite('non_jumping_beads_%04d.tif' % frame_index, is_non_jumping_bead_image)
|
|
# extract the surface points
|
|
surface_y = np.ndarray(shape=(is_non_jumping_bead_image.shape[1],), dtype=int)
|
|
surface_y.fill(is_non_jumping_bead_image.shape[0])
|
|
contours, hierarchy = cv2.findContours(is_non_jumping_bead_image, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
|
|
print('num contours: ', len(contours))
|
|
assert(len(contours) == 1)
|
|
contour = contours[0]
|
|
for point in contour:
|
|
x, y = point[0]
|
|
surface_y[x] = min(surface_y[x], y)
|
|
print(surface_y)
|
|
surface_pts = []
|
|
for x in range(len(surface_y)):
|
|
y = surface_y[x]
|
|
if y != is_non_jumping_bead_image.shape[0]:
|
|
surface_pts.append((x, y))
|
|
print(surface_pts)
|
|
x = [pt[0] for pt in surface_pts]
|
|
y = [pt[1] for pt in surface_pts]
|
|
lin_regress_result = scipy.stats.linregress(x, y)
|
|
print(lin_regress_result)
|
|
|
|
@staticmethod
|
|
def create_isbead_image(trac_data: TracData, frame_index: int, particle_radius: float):
|
|
image = np.zeros(shape=trac_data.image_size, dtype=float)
|
|
print(trac_data.pts)
|
|
frame_pts = trac_data.pts[trac_data.pts[:, 0] == frame_index]
|
|
for x in frame_pts[:, 2:4]:
|
|
center_coordinates = (int(x[0]), int(x[1]))
|
|
# center_coordinates = (10,10)
|
|
print(center_coordinates)
|
|
cv2.circle(image, center_coordinates, int(particle_radius), color=1.0, thickness=-1)
|
|
return image
|
|
|
|
|
|
def main():
|
|
# python3 ./tractrac.git/Python/tractrac.py -f ./grassloper.git/samples/sample002.avi --output 1 -a --saveplot
|
|
trac_data = hdf5_to_trac_data('./grassloper.git/samples/TracTrac/sample002_track.hdf5')
|
|
slope_finder = SlopeFinder(beads_radius=11.0)
|
|
for frame_index in range(1, 2):
|
|
slope_finder.find_slope(trac_data, frame_index=frame_index)
|
|
|
|
# slope_image = create_slope_image(image_width=128, image_height=256, p1_angle=0.7, p1_radius=3.0)
|
|
# cv2.imwrite(str(Path('toto.tif')), slope_image)
|
|
|
|
|
|
if __name__ == '__main__':
|
|
main()
|